SCM

SCM Repository

[returnanalytics] View of /pkg/Meucci/R/EntropyProg.R
ViewVC logotype

View of /pkg/Meucci/R/EntropyProg.R

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3717 - (download) (annotate)
Tue Jun 23 08:26:43 2015 UTC (2 years, 1 month ago) by xavierv
File size: 14325 byte(s)
- fixed HermiteGrid_demo demo script
#' @title Compute posterior (=change of measure) with Entropy Pooling.
#'
#' @description
#' Entropy Pooling program for blending views on scenarios with a 
#' prior scenario-probability distribution.
#'
#' @details
#' The entropy program will change the initial predictive distribution 'p' to a
#' new set 'p_' that satisfies specified moment conditions but changes other
#' properties of the new distribution, the least by minimizing the relative
#' entropy between the two distributions. Theoretical note: Relative Entropy
#' (Kullback-Leibler information criterion KLIC) is an asymmetric measure. 
#'
#' We retrieve a new set of probabilities for the joint-scenarios using the
#' Entropy pooling method of the many choices of 'p' that satisfy the views,
#' we choose a 'p' that minimize the entropy or distance of the new probability
#' distribution to the prior joint-scenario probabilities, we use the
#' Kullback-Leibler divergence or the relative entropy dist(p,q): Sum across all
#' scenarios [p-t * ln( p-t / q-t)]. Therefore we define solution as
#' p* = argmin (choice of p) [sum across all scenarios: p-t * ln( p-t / q-t)],
#' such that 'p' satisfies views. The views modify the prior in a coherent
#' manner (minimizing distortion). We formulate the stress tests of the baseline
#' scenarios as linear constraints on yet-to-be defined probabilities.
#'
#' Note that the numerical optimization acts on a very limited number of
#' variables equal to the number of views. It does not act directly on the very
#' large number of variables of interest, namely the probabilities of the Monte
#' Carlo scenarios. This feature guarantees the numerical feasability of entropy
#' optimization. 
#'
#' Note that new probabilities are generated in much the same way that the 
#' state-price density modifies objective probabilities of pay-offs to
#' risk-neutral probabilities in contingent-claims asset pricing.
#'
#'
#' @param  p        a vector of initial probabilities based on prior (reference
#'                  model, empirical distribution, etc.). Sum of 'p' must be 1
#' @param  Aeq      matrix consisting of equality constraints (paired with
#'                  argument 'beq'). Denoted as 'H' in the Meucci paper. 
#'                  (denoted as 'H' in the "Meucci - Flexible Views Theory & 
#'                  Practice" paper formlua 86 on page 22)
#' 
#' @param  beq      vector corresponding to the matrix of equality constraints
#'                  (paired with argument 'Aeq'). Denoted as 'h' in the paper
#' @param  A        matrix consisting of inequality constraints (paired with
#'                  argument 'b'). Denoted as 'F' in the Meucci paper
#' @param  b        vector consisting of inequality constraints (paired with
#'                  matrix A). Denoted as 'f' in the Meucci paper
#'
#' \deqn{ \tilde{p}  \equiv  argmin_{Fx \leq f, Hx  \equiv  h}  \big\{ \sum_1^J
#'  x_{j}  \big(ln \big( x_{j} \big) - ln \big( p_{j} \big) \big)  \big\} \\ 
#' \ell  \big(x,  \lambda,  \nu \big)  \equiv  x'  \big(ln \big(x\big) - ln 
#' \big(p\big) \big) +   \lambda' \big(Fx - f\big)  +   \nu' \big(Hx - h\big)}
#'
#' @return   p_                        revised probabilities based on entropy
#'                                     pooling
#' @return   optimizationPerformance   a list with status of optimization,
#'                                     value, number of iterations and sum of
#'                                     probabilities.
#'
#' @author Ram Ahluwalia \email{ram@@wingedfootcapital.com}
#'
#' @references 
#' A. Meucci - "Fully Flexible Views: Theory and Practice". See page 22 for 
#' illustration of numerical implementation
#'
#' Symmys site containing original MATLAB code \url{http://www.symmys.com}
#'
#' NLOPT open-source optimization site containing background on algorithms 
#' \url{http://ab-initio.mit.edu/wiki/index.php/NLopt}
#'
#' We use the information-theoretic estimator of Kitamur and Stutzer (1997). 
#'
#' Reversing 'p' and 'p_' leads to the empirical likelihood" estimator of Qin 
#' and Lawless (1994). 
#'
#' See Robertson et al, "Forecasting Using Relative Entropy" (2002) for theory
#'
#' @export

EntropyProg <- function(p, A = NULL, b = NULL, Aeq, beq) {

  if(!length(A))
    A <- matrix(,nrow = 0, ncol = 0)
  if(!length(b))
    b <- matrix(,nrow = 0, ncol = 0)

  # count the number of constraints
  # K_ is the number of inequality constraints in the matrix-vector pair A-b
  K_ <- nrow(A)
  # K is the number of equality views in the matrix-vector pair Aeq-beq
  K  <- nrow(Aeq)

  # parameter checks
  if (K_ + K == 0)
    stop("at least one equality or inequality constraint must be specified")

  if (((.999999 < sum(p)) & (sum(p) < 1.00001)) == FALSE)
    stop("sum of probabilities from prior distribution must equal 1")

  if (nrow(Aeq) != length(beq)){
    stop("number of inequality constraints in matrix Aeq must match
      number of elements in vector beq")
  }
  if (nrow(A) != length(b)) {
    stop("number of equality constraints in matrix A must match number of
          elements in vector b")
  }

  # calculate derivatives of constraint matrices
  A_   <- t(A)
  Aeq_ <- t(Aeq)
  beq_ <- t(beq)

  # starting guess for optimization search with length = to number of views
  x0 <- matrix(0, nrow = K_ + K, ncol = 1)

  # equality constraints only
  if (!K_) {
    # define maximum likelihood, gradient, and hessian functions for
    # unconstrained and constrained optimization
    eval_f_list <- function(v) {
      # cost function for unconstrained optimization
      # (no inequality constraints)
      # L:     Lagrangian dual function (without inequality constraints)
      #        See formula 88 on p. 22 in "Meucci - Fully Flexible Views -
      #        Theory and Practice (2010)".
      # t(x):  the derivative x'
      # Aeq_:  the derivative of the Aeq matrix of equality constraints
      #        (denoted as 'H in the paper)
      # beq_:  transpose of the vector associated with Aeq equality
      #        constraints
      # L=  x'  *  (log(x) - log(p) + Aeq_  *  v) -   beq_ *  v
      #    1xJ  *   (Jx1   - Jx1  + JxN+1 * N+1x1) - 1xN+1 * N+1x1

      x <- exp(log(p) - 1 - Aeq_ %*% v)
      x <- apply(cbind(x, 10 ^ -32), 1, max)  # robustification
      L <- t(x) %*% (log(x) - log(p) + Aeq_ %*% v) - beq_ %*% v
      mL <- -L # take negative values since we want to maximize
      # evaluate gradient
      gradient <- beq - Aeq %*% x

      # evaluate Hessian
      # We comment this out for now -- to be used when we find an
      # optimizer that can utilize the Hessian in addition to the gradient
      # H = (Aeq %*% ((x %*% ones(1,K)) * Aeq_))
      # Hessian computed by Chen Qing, Lin Daimin, Meng Yanyan,Wang Weijun

      return(list(objective = mL, gradient = gradient))
    }

    # setup unconstrained optimization
    start <- Sys.time()
    opts <- list(algorithm = "NLOPT_LD_LBFGS", xtol_rel = 1.0e-6,
                check_derivatives = TRUE, check_derivatives_print = "all",
                print_level = 2, maxeval = 1000)
    optimResult <- nloptr(x0 = x0, eval_f = eval_f_list, opts = opts)
    end <- Sys.time()
    print(c("Optimization completed in ", end - start))
    rm(start)
    rm(end)

    if (optimResult$status < 0) {
      print(c("Exit code ", optimResult$status))
      stop("Error: The optimizer did not converge")
    }

    # return results of optimization
    v <- optimResult$solution
    p_ <- exp(log(p) - 1 - Aeq_ %*% v)
    optimizationPerformance <- list(converged = (optimResult$status > 0),
                                    ml = optimResult$objective,
                                    iterations = optimResult$iterations,
                                    sumOfProbabilities = sum(p_))
  } else {
    #--- case inequality constraints are specified

    #-- setup variables for constrained optimization

    # -1 * Identity Matrix with dimension equal to number of constraints
    InqMat <- -diag(1, K_ + K)
    # drop rows corresponding to equality constraints
    InqMat <- InqMat[-c(K_ + 1:nrow(InqMat)),]
    InqVec <- matrix(0, K_, 1)

    #-- define maximum likelihood, gradient, and hessian functions for
    # constrained optimization.

    # function used to evalute A %*% x <= 0 or A %*% x <= c(0,0) if there is
    # more than one inequality constraint
    InqConstraint <- function(x) {
      return(InqMat %*% x)
    }

    jacobian_constraint <- function(x) {
      return(InqMat)
    }
    # Jacobian of the constraints matrix. One row per constraint, one column
    # per control parameter (x1,x2)
    # Turns out the Jacobian of the constraints matrix is always equal to
    # InqMat

    nestedfunC <- function(lv) {
      lv <- as.matrix(lv)
      # inequality Lagrange multiplier
      l <- lv[1:K_, , drop = FALSE]
      # equality lagrange multiplier
      v <- lv[(K_ + 1):length(lv), , drop = FALSE]
      x <- exp(log(p) - 1 - A_ %*% l - Aeq_ %*% v)
      x <- apply(cbind(x, 10 ^ -32), 1, max)

      # L is the cost function used for constrained optimization.
      # L is the Lagrangian dual function with inequality constraints and
      # equality constraints
      L <- t(x) %*% (log(x) - log(p)) + t(l) %*% (A %*% x - b) + t(v) %*%
           (Aeq %*% x - beq)
      objective <- -L  # take negative values since we want to maximize

      # calculate the gradient
      gradient <- rbind(b - A %*% x, beq - Aeq %*% x)

      ##-- compute the Hessian (commented out since no R optimizer
      ## supports use of Hessian)
      ## Hessian by Chen Qing, Lin Daimin, Meng Yanyan, Wang Weijun
      # onesToK_ = array(rep(1, K_))
      # onesToK = array(rep(1, K))
      # x = as.matrix(x)
      # H11 = A %*% ((x %*% onesToK_) * A_)
      # H12 = A %*% ((x %*% onesToK) * Aeq_)
      # H21 = Aeq %*% ((x %*% onesToK_) * A_)
      # H22 = Aeq %*% ((x %*% onesToK) * Aeq_)
      # H1 = cbind(H11, H12)
      # H2 = cbind(H21, H22)
      ## H: Hessian for constrained optimization
      # H = rbind(H1, H2)

      return(list(objective = objective, gradient = gradient))
    }

    # find minimum of constrained multivariate function
    start <- Sys.time()
    # Note: other candidates for constrained optimization in library nloptr:
    # NLOPT_LD_SLSQP, NLOPT_LD_MMA, NLOPT_LN_AUGLAG, NLOPT_LD_AUGLAG_EQ
    # See NLOPT open-source site for more details:
    # http://ab-initio.mit.edu/wiki/index.php/NLopt
    local_opts <- list(algorithm = "NLOPT_LD_SLSQP", xtol_rel = 1.0e-6,
                       check_derivatives = TRUE, 
                       check_derivatives_print = "all", eval_f = nestedfunC,
                       eval_g_ineq = InqConstraint,
                       eval_jac_g_ineq = jacobian_constraint)
    optimResult <- nloptr(x0 = x0, eval_f = nestedfunC,
                          eval_g_ineq = InqConstraint,
                          eval_jac_g_ineq = jacobian_constraint,
                          opts = list(algorithm = "NLOPT_LD_AUGLAG",
                                      local_opts = local_opts,
                                      print_level = 2, maxeval = 1000,
                                      check_derivatives = TRUE,
                                      check_derivatives_print = "all",
                                      xtol_rel = 1.0e-6))
    end <- Sys.time()
    print(c("Optimization completed in ", end - start))
    rm(start)
    rm(end)

    if (optimResult$status < 0) {
      print(c("Exit code ", optimResult$status))
      stop("Error: The optimizer did not converge")
    }

    # return results of optimization

    lv <- matrix(optimResult$solution, ncol = 1)
    # inequality Lagrange multipliers
    l <- lv[1:K_, , drop = FALSE]
    # equality Lagrange multipliers
    v <- lv[(K_ + 1):nrow(lv), , drop = FALSE]
    p_ <- exp(log(p) - 1 - A_ %*% l - Aeq_ %*% v)
    optimizationPerformance <- list(converged = (optimResult$status > 0),
                                    ml = optimResult$objective,
                                    iterations = optimResult$iterations,
                                    sumOfProbabilities = sum(p_))
  }

  print(optimizationPerformance)

  if (sum(p_) < .999)
    stop("Sum or revised probabilities is less than 1!")
  if (sum(p_) > 1.001)
    stop("Sum or revised probabilities is greater than 1!")

  return (list (p_ = p_, optimizationPerformance = optimizationPerformance))
}



#' Generates histogram
#'
#' @param X       vector containing the data points
#' @param p       vector containing the probabilities for each of the data 
#'                points in X
#' @param nBins   expected number of Bins the data set is to be broken down into
#' @param freq    boolean variable to indicate whether the graphic is a 
#'                representation of frequencies
#' @param main    title for the plot
#' @param xlim    limits for the x-axis
#' @param ylim    limits for the y-axis
#'
#' @return a list with 
#'             f   the frequency for each midpoint
#'             x   the midpoints of the nBins intervals
#'
#' @references 
#' \url{http://www.symmys.com}
#' See Meucci script pHist.m used for plotting
#' @author Ram Ahluwalia \email{ram@@wingedfootcapital.com} and Xavier Valls 
#' \email{xaviervallspla@@gmail.com}
#'
#' @export

PHist <- function(X, p, nBins, freq = FALSE,
                  main = "Portfolio return distribution",
                  xlim = NULL, ylim = NULL) {

  if (length(match.call()) < 3) {
    J <- dim(X)[1]
    nBins <- round(10 * log(J))
  }

  x <- hist(x = X, breaks = nBins, plot = FALSE)$breaks
  D <- x[2] - x[1]

  N <- length(x)
  np <- zeros(N, 1)

  for (s in 1:N) {
    # The boolean Index is true is X is within the interval centered at x(s) and
    # within a half-break distance
    Index <- (X >= x[s] - D / 2) & (X <= x[s] + D / 2)
    # np = new probabilities?
    np[s] <- sum(p[Index])
    f <- np / D
  }

  plot(x, f, type = "h", main = main, xlim = xlim, ylim = ylim)

  return(list(f = f, x = x))
}

R-Forge@R-project.org
ViewVC Help
Powered by ViewVC 1.0.0  
Thanks to:
Vienna University of Economics and Business University of Wisconsin - Madison Powered By FusionForge