[Returnanalytics-commits] r3895 - in pkg/Meucci: . R data demo man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sun Aug 2 19:21:04 CEST 2015
Author: xavierv
Date: 2015-08-02 19:21:04 +0200 (Sun, 02 Aug 2015)
New Revision: 3895
Added:
pkg/Meucci/R/DynamicPortfolioManagement.R
pkg/Meucci/R/MVOUPosterior.R
pkg/Meucci/R/MVOUPrior.R
pkg/Meucci/R/QuadraticMatVC.R
pkg/Meucci/data/dynamicManagement.rda
pkg/Meucci/demo/S_DynamicManagementCase1.R
pkg/Meucci/demo/S_DynamicManagementCase2.R
pkg/Meucci/man/BellmanEq_CS1.Rd
pkg/Meucci/man/BellmanEq_CS2.Rd
pkg/Meucci/man/MVOU_Posterior.Rd
pkg/Meucci/man/MVOU_Prior.Rd
pkg/Meucci/man/QuadraticMat_Vc.Rd
pkg/Meucci/man/dynamicManagement.Rd
Modified:
pkg/Meucci/DESCRIPTION
pkg/Meucci/NAMESPACE
pkg/Meucci/R/FullyFlexibleBayesNets.R
pkg/Meucci/R/data.R
pkg/Meucci/demo/FullyFlexibleBayesNets.R
pkg/Meucci/man/Equities.Rd
pkg/Meucci/man/JGB.Rd
pkg/Meucci/man/StockSeries.Rd
pkg/Meucci/man/TimeSeries.Rd
pkg/Meucci/man/UsSwapRates.Rd
pkg/Meucci/man/bondAttribution.Rd
pkg/Meucci/man/butterfliesAnalytics.Rd
pkg/Meucci/man/covNRets.Rd
pkg/Meucci/man/db.Rd
pkg/Meucci/man/dbFFP.Rd
pkg/Meucci/man/db_FX.Rd
pkg/Meucci/man/derivatives.Rd
pkg/Meucci/man/fILMR.Rd
pkg/Meucci/man/factorsDistribution.Rd
pkg/Meucci/man/fixedIncome.Rd
pkg/Meucci/man/freaqEst.Rd
pkg/Meucci/man/highYieldIndices.Rd
pkg/Meucci/man/implVol.Rd
pkg/Meucci/man/linRet.Rd
pkg/Meucci/man/linearModel.Rd
pkg/Meucci/man/returnsDistribution.Rd
pkg/Meucci/man/sectorsSnP500.Rd
pkg/Meucci/man/sectorsTS.Rd
pkg/Meucci/man/securitiesIndustryClassification.Rd
pkg/Meucci/man/securitiesTS.Rd
pkg/Meucci/man/swap2y4y.Rd
pkg/Meucci/man/swapParRates.Rd
pkg/Meucci/man/swaps.Rd
Log:
Added Dynamic Portfolio Management paper scripts and case studies, reformatted and modified FullyFlexibleBayesNets scripts
Modified: pkg/Meucci/DESCRIPTION
===================================================================
--- pkg/Meucci/DESCRIPTION 2015-07-31 09:27:25 UTC (rev 3894)
+++ pkg/Meucci/DESCRIPTION 2015-08-02 17:21:04 UTC (rev 3895)
@@ -29,7 +29,6 @@
R (>= 2.14.0),
zoo,
xts (>= 0.8),
- matlab,
pracma,
R.utils,
mvtnorm,
@@ -38,6 +37,7 @@
kernlab,
nloptr,
limSolve,
+ linprog,
Suggests:
Matrix,
MASS,
Modified: pkg/Meucci/NAMESPACE
===================================================================
--- pkg/Meucci/NAMESPACE 2015-07-31 09:27:25 UTC (rev 3894)
+++ pkg/Meucci/NAMESPACE 2015-08-02 17:21:04 UTC (rev 3895)
@@ -1,5 +1,7 @@
# Generated by roxygen2 (4.1.1): do not edit by hand
+export(BellmanEq_CS1)
+export(BellmanEq_CS2)
export(BlackLittermanFormula)
export(BlackScholesCallPrice)
export(BlackScholesCallPutPrice)
@@ -37,6 +39,8 @@
export(LognormalMoments2Parameters)
export(LognormalParam2Statistics)
export(LongShortMeanCVaRFrontier)
+export(MVOU_Posterior)
+export(MVOU_Prior)
export(MaxRsqCS)
export(MaxRsqTS)
export(MleRecursionForStudentT)
@@ -55,6 +59,7 @@
export(PlotVolVsCompositionEfficientFrontier)
export(Prior2Posterior)
export(ProjectionStudentT)
+export(QuadraticMat_Vc)
export(QuantileMixture)
export(RIEfficientFrontier)
export(RandNormalInverseWishart)
Added: pkg/Meucci/R/DynamicPortfolioManagement.R
===================================================================
--- pkg/Meucci/R/DynamicPortfolioManagement.R (rev 0)
+++ pkg/Meucci/R/DynamicPortfolioManagement.R 2015-08-02 17:21:04 UTC (rev 3895)
@@ -0,0 +1,384 @@
+#' Solves the Bellman Equation for the case study 1.
+#'
+#' @details In Case Study 1 there is only one risk driver (the 10 year rate) and
+#' only one view. The view is that the expected value of the 10-year rate will
+#' be the actual value minus 50 basis points at t^view = 1 year from the current
+#' time.
+#' The solution is analytical in the prior case.
+#' The solution is found recursivelly in the posterior case starting from
+#' t_view and going back to time 0 with a time step = tau.
+#'
+#' In case study 1: n_ = 1 N_meanViews = 1
+#'
+#' @param eta [scalar] overall weight of the market impact of transactions
+#' @param gamma [scalar] risk aversion parameter
+#' @param lambda [scalar] discounting parameter
+#' @param tau [scalar] trading interval
+#' @param theta [n_ x n_] transition matrix of the MVOU process
+#' @param mu [n_ x 1] drift vector of the MVOU process
+#' @param sig2 [n_ x n_] covariance parameters of the MVOU process
+#' @param c2 [n_ x n_] matrix of the market impact
+#' @param b_legacy [n_ x 1] legacy portfolio exposure at time 0
+#' @param x [t_ x n_] path of the risk drivers (with time step = tau)
+#' @param t_view [1 x N_MeanViews] times of the views
+#' @param view [1 x N_MeanViews] views on the risk drivers
+#'
+#' @return prior [t_x n_ matrix] optimal prior exposure
+#' @return post [t_x n_ matrix] optimal posterior exposure
+#'
+#' @references
+#' A. Meucci - "Dynamic Portfolio Management with Views at Multiple Horizons"
+#' \url{http://symmys.com/node/831}. See Meucci script for "BellmanEq_CS1.m"
+#'
+#' @author Xavier Valls \email{xavievallspla@@gmail.com}
+#' @export
+
+BellmanEq_CS1 <- function(eta, gamma, lambda, tau, theta, mu, sig2, c2,
+ b_legacy, x, t_view, view) {
+ t_ <- nrow(x)
+ n_ <- length(theta)
+
+ ##############################################################################
+ #compute the prior at time 0
+ Prior0 <- MVOU_Prior(c(0, tau), x[1], theta, sig2, mu)
+ # first period covariance matrix
+ sig2_1 <- Prior0$cov[(n_ + 1):(2 * n_), (n_ + 1):(2 * n_)]
+
+ ##############################################################################
+ #Coefficients of the Bellman equation according to the prior
+
+ alpha_prior <- Prior0$mean_cost[(n_ + 1):(2 * n_)]
+ beta_prior <- Prior0$mean_lin[(n_ + 1):(2 * n_), 1:n_] - diag(1, n_)
+ HATsig2 <- exp(lambda) * (eta * c2) ^ (-1 / 2) * gamma * sig2_1 * (eta * c2) ^
+ (-1 / 2)
+ HATpsi_bb <- (0.25 * ( HATsig2 + diag(1, n_) * (exp(lambda) - 1)) ^ 2 +
+ HATsig2) ^ (1 / 2) - 0.5 * (HATsig2 + diag(1, n_) *
+ (exp(lambda) - 1))
+ psi_bb_prior <- (eta * c2) ^ (1 / 2) * HATpsi_bb * (eta * c2) ^ (1 / 2)
+ q_prior <- gamma * sig2_1 + eta * c2 + exp(-lambda) * psi_bb_prior
+ tmp <- (eta * c2 * (solve(q_prior) %*% beta_prior))
+ psi_bx_prior <- solve(diag(1, n_ ^ 2) - exp(-lambda) * kron(t(beta_prior) +
+ diag(1, n_), eta * ( c2 / q_prior))) %*% array(tmp)
+ psi_bx_prior <- matrix(psi_bx_prior, nrow = n_)
+ psi_b_prior <- solve(q_prior / (eta * c2) - exp(-lambda) * diag(1, n_)) %*%
+ (diag(1, n_) + exp(-lambda) * psi_bx_prior) * alpha_prior
+
+ # #Alternatively, if and only if c2 <- sig2_1
+ # a_ <- (sqrt(4*gamma*eta*exp(-lambda)+(gamma+(1-exp(-lambda))*eta)^2) -
+ # ((1-exp(-lambda))*eta+gamma))/(2*exp(-lambda))
+ # psi_bb_prior <- a_*sig2_1
+ # psi_bx_prior <- eta*(beta_prior)*inv((gamma+eta+exp(-lambda)*a_)*
+ # diag(1, n_)-eta*exp(-lambda)*(beta_prior+diag(1, n_)))
+ # psi_b_prior <- eta*(alpha_prior + exp(-lambda)*psi_bx_prior*alpha_prior)/
+ # (gamma+eta+exp(-lambda)*a_-exp(-lambda)*eta)
+
+ ##############################################################################
+ #Coefficients of the Bellman equation according to the posterior distribution
+ ##############################################################################
+
+ ##############################################################################
+ #Inizialize the variables
+ Hor <- ceil(max(t_view) / tau)
+ psi_t_bb <- array(0, dim = c(n_,n_, Hor))
+ q_t <- array(0, dim = c(n_,n_, Hor))
+ psi_t_bx <- array(0, dim = c(n_,n_, Hor))
+ psi_t_b <- matrix(0, n_, Hor)
+ alpha_t <- matrix(0, n_, Hor)
+ beta_t <- array(0, dim = c(n_,n_, Hor))
+
+ ##############################################################################
+ #Set the boundary conditions asintotically. After the last view, the
+ #solution is equal to the prior
+ psi_t_bb[1:n_, 1:n_, Hor] <- psi_bb_prior
+ q_t[1:n_, 1:n_, Hor] <- q_prior
+ psi_t_bx[1:n_, 1:n_, Hor] <- psi_bx_prior
+ psi_t_b[1:n_, Hor] <- psi_b_prior
+ alpha_t[1:n_, Hor] <- alpha_prior
+ beta_t[1:n_, 1:n_, Hor] <- beta_prior
+
+ for (k in seq(Hor - 1, 1, -1)) {
+ q <- drop(q_t[1:n_, 1:n_, k + 1])
+ alpha <- drop(alpha_t[1:n_, k + 1])
+ beta <- drop(beta_t[1:n_, 1:n_, k + 1])
+ psi_bx <- drop(psi_t_bx[1:n_, 1:n_, k + 1])
+ psi_b <- drop(psi_t_b[1:n_, k + 1])
+
+ ############################################################################
+ #update of the coefficients of the Bellman equation
+ psi_t_b[1:n_, k] <- eta * c2 * (solve(q) %*% (alpha + exp(-lambda) *
+ psi_bx * alpha + exp(-lambda) * psi_b))
+ psi_t_bx[1:n_, 1:n_, k] <- eta * c2 * (solve(q) %*% (beta + exp(-lambda) *
+ psi_bx * (beta + diag(1, n_))))
+ psi_t_bb[1:n_, 1:n_, k] <- -eta ^ 2 * c2 * (solve(q) %*% c2) + eta * c2
+
+ ############################################################################
+ #set the monitoring times of interest
+ t <- c(0, tau, t_view - k * tau)
+ if ((t[length(t)] - t[length(t) - 1]) < t[length(t)] * 10 ^ (-10))
+ t <- t[-length(t)]
+ T_ <- length(t)
+
+ #compute the prior
+ Prior <- MVOU_Prior(t, 0, theta, sig2, mu)
+
+ #set the views
+ grid <- meshgrid(t, 1:n_)
+ T <- grid$X
+ N <- grid$Y
+ N_Meanviews <- 1 #Number of views on expectations
+ v_mu_tmp <- array(0, dim = c(N_Meanviews, n_, T_))
+ v_mu_tmp[1, 1, T_] <- 1
+ mu_view <- view[1]
+ v_mu <- matrix(v_mu_tmp, nrow = nrow(v_mu_tmp))
+ views <- list()
+ views$N_Meanviews <- N_Meanviews
+ views$N_Covviews <- c()
+ views$dimension <- array(N)
+ views$monitoring_time <- array(T)
+ views$v_mu <- v_mu
+ views$v_sig <- NaN
+ views$mu_view <- mu_view
+ views$sig2_view <- c()
+
+ ############################################################################
+ #update the Posterior moments of the process
+ Posterior <- MVOU_Posterior(Prior, views)
+ alpha_t[1:n_, k] <- Posterior$mean_cost[(n_ + 1):(2 * n_)]
+ beta_t[1:n_, 1:n_, k] <- Posterior$mean_lin[(n_ + 1):(2 * n_), 1:n_] -
+ diag(1, n_)
+ sig2_1 <- Posterior$cov[(n_ + 1):(2 * n_), (n_ + 1):(2 * n_)]
+
+ ############################################################################
+ #update matrix q_t
+ psi_bb <- drop(psi_t_bb[1:n_, 1:n_, k])
+ q_t[1:n_, 1:n_, k] <- gamma * sig2_1 + eta * c2 + exp(-lambda) * psi_bb
+ }
+
+ b_prior <- matrix(NaN, t_, 1 )
+ b_post <- matrix(NaN, t_, 1 )
+
+ #Reconstructing the optimal exposure on the simulated path x
+ for (t in 1:t_) {
+ if (t == 1) {
+ b_legacy_prior <- b_legacy
+ b_legacy_post <- b_legacy
+ } else {
+ b_legacy_prior <- b_prior[t - 1]
+ b_legacy_post <- b_post[t - 1]
+ }
+ # prior exposure
+ b_prior[t] <- solve(q_prior) %*% (alpha_prior + (beta_prior) * x[t] + eta *
+ c2 * b_legacy_prior + exp(-lambda) * psi_bx_prior * (alpha_prior +
+ (beta_prior + diag(1, n_)) * x[t]) + exp(-lambda) * psi_b_prior)
+ #posterior exposure
+ q <- drop(q_t[1:n_, 1:n_, t])
+ alpha <- drop(alpha_t[1:n_, t])
+ beta <- drop(beta_t[1:n_, 1:n_, t])
+ psi_bx <- drop(psi_t_bx[1:n_, 1:n_, t])
+ psi_b <- drop(psi_t_b[1:n_, t])
+ l <- alpha + beta * x[t] + eta * c2 * b_legacy_post + exp(-lambda) *
+ psi_bx * (alpha + (beta + diag(1, n_)) * x[t]) + exp(-lambda) * psi_b
+ b_post[t] <- solve(q) %*% l
+ }
+ return(list(prior = b_prior, post = b_post))
+}
+
+#' Solves the Bellman Equation for the case study 2.
+#'
+#' @details In Case Study 2 we consider two risk drivers, the 10 year rate and
+#' the TIP spread, and two non-synchronous views on them. The view on the rate
+#' is that its expected value will be the actual value minus 50 basis points
+#' at t_viewX = 1 year from the current time (as in Case Study 1).
+#' The view on the TIP spread is that its expected value will be the actual
+#' value plus 50 basis points at t_viewTIP = 0.75 years.
+#'
+#' In case study 2: n_ = 2 k_ = 1 N_meanViews = 2
+#'
+#' @param eta [scalar] overall weight of the market impact of transactions
+#' @param gamma [scalar] risk aversion parameter
+#' @param lambda [scalar] discounting parameter
+#' @param tau [scalar] trading interval
+#' @param theta [n_ x n_] transition matrix of the MVOU process
+#' @param mu [n_ x 1] drift vector of the MVOU process
+#' @param sig2 [n_ x n_] covariance parameters of the MVOU process
+#' @param c2 [n_ x n_] matrix of the market impact
+#' @param b_legacy [n_ x 1] legacy portfolio exposure at time 0
+#' @param x [t_ x n_] path of the risk drivers (with time step = tau)
+#' @param t_view [1 x N_MeanViews] times of the views
+#' @param view [1 x N_MeanViews] views on the risk drivers
+#' @param i_view [1 x N_MeanViews] vector of the labels of the risk drivers
+#' to which views refer
+#' @param omega [k_ x n_] matrix to select the investible risk drivers
+#'
+#' @return prior [t_x n_ matrix] optimal prior exposure
+#' @return post [t_x n_ matrix] optimal posterior exposure
+#'
+#' @references
+#' A. Meucci - "Dynamic Portfolio Management with Views at Multiple Horizons"
+#' \url{http://symmys.com/node/831}. See Meucci script for "BellmanEq_CS2.m"
+#'
+#' @author Xavier Valls \email{xavievallspla@@gmail.com}
+#' @export
+
+BellmanEq_CS2 <- function(eta, gamma, lambda, tau, theta, mu, sig2, c2,
+ b_legacy, x, t_view, view, i_view, omega) {
+
+ t_ <- nrow(x)
+ n_ <- nrow(theta)
+ k_ <- length(c2)
+
+ ##############################################################################
+ #compute the prior at time 0
+ Prior0 <- MVOU_Prior(c(0, tau), matrix(x[1:n_]), theta, sig2, mu)
+ # first period covariance matrix
+ sig2_1 <- Prior0$cov[(n_ + 1):(2 * n_), (n_ + 1):(2 * n_)]
+ print(omega)
+ print(sig2_1)
+ print(t(omega))
+ sig2_1 <- omega %*% sig2_1 %*% t(omega)
+ ##############################################################################
+ #Coefficients of the Bellman equation according to the prior
+
+ alpha_prior <- Prior0$mean_cost[(n_ + 1):(2 * n_)]
+ beta_prior <- Prior0$mean_lin[(n_ + 1):(2 * n_), 1:n_] - diag(1, n_)
+ HATsig2 <- exp(lambda) * (eta * c2) ^ (-1 / 2) * gamma * sig2_1 * (eta * c2) ^
+ (-1 / 2)
+ HATpsi_bb <- (0.25 * ( HATsig2 + diag(1, k_) * (exp(lambda) - 1)) ^ 2 +
+ HATsig2) ^ (1 / 2) - 0.5 * (HATsig2 + diag(1, k_) *
+ (exp(lambda) - 1))
+ psi_bb_prior <- (eta * c2) ^ (1 / 2) * HATpsi_bb * (eta * c2) ^ (1 / 2)
+ q_prior <- gamma * sig2_1 + eta * c2 + exp(-lambda) * psi_bb_prior
+ tmp <- (eta * c2 * (solve(q_prior) %*% (omega %*% beta_prior)))
+ psi_bx_prior <- solve(diag(1, k_ * n_) - exp(-lambda) * kron(t(beta_prior) +
+ diag(1, n_), eta * ( c2 / q_prior))) %*% array(tmp)
+ psi_bx_prior <- matrix(psi_bx_prior, nrow = k_, ncol = n_)
+ psi_b_prior <- solve(q_prior / (eta * c2) - exp(-lambda) * diag(1, k_)) %*%
+ (omega + exp(-lambda) * psi_bx_prior) %*% alpha_prior
+
+ ##############################################################################
+ #Coefficients of the Bellman equation according to the posterior distribution
+ ##############################################################################
+
+ ##############################################################################
+ #Inizialize the variables
+ Hor <- ceil(max(t_view) / tau)
+ psi_t_bb <- array(0, dim = c( k_, k_, Hor))
+ q_t <- array(0, dim = c( k_, k_, Hor))
+ psi_t_bx <- array(0, dim = c( k_, n_, Hor))
+ psi_t_b <- matrix(0, k_, Hor)
+ alpha_t <- matrix(0, n_, Hor)
+ beta_t <- array(0, dim = c( n_, n_, Hor))
+
+ ##############################################################################
+ #Set the boundary conditions asintotically. After the last view, the
+ #solution is equal to the prior
+ psi_t_bb[1:k_, 1:k_, Hor] <- psi_bb_prior
+ q_t[1:k_, 1:k_, Hor] <- q_prior
+ psi_t_bx[1:k_, 1:n_, Hor] <- psi_bx_prior
+ psi_t_b[1:k_, Hor] <- psi_b_prior
+ alpha_t[1:n_, Hor] <- alpha_prior
+ beta_t[1:n_, 1:n_, Hor] <- beta_prior
+
+ for (k in seq(Hor - 1, 1, -1)) {
+ q <- drop(q_t[1:k_, 1:k_, k + 1])
+ alpha <- drop(alpha_t[1:n_, k + 1])
+ beta <- drop(beta_t[1:n_, 1:n_, k + 1])
+ psi_bx <- drop(psi_t_bx[1:k_, 1:n_, k + 1])
+ psi_b <- drop(psi_t_b[1:k_, k + 1])
+
+ ############################################################################
+ #update of the coefficients of the Bellman equation
+ psi_t_b[1:k_, k] <- eta * c2 * (solve(q) %*% (omega %*% alpha +
+ exp(-lambda) * psi_bx %*% alpha + exp(-lambda) * psi_b))
+ psi_t_bx[1:k_, 1:n_, k] <- eta * c2 * (solve(q) %*% (omega %*% beta +
+ exp(-lambda) * psi_bx %*% (beta + diag(1, n_))))
+ psi_t_bb[1:k_, 1:k_, k] <- -eta ^ 2 * c2 * (solve(q) %*% c2) + eta * c2
+
+ ############################################################################
+ #set the monitoring times of interest
+ t <- c(0, tau, t_view - k * tau)
+ t <- sort(t)
+ t <- t[t>=0]
+ idx <- which(diff(t) < tau * 10 ^ -10)
+ t <- t(setdiff(1:length(t), idx))
+ T_ <- length(t)
+
+ #compute the prior
+ Prior <- MVOU_Prior(t, matrix(c(0,0)), theta, sig2, mu)
+
+ #set the views
+ grid <- meshgrid(t, 1:n_)
+ T <- grid$X
+ N <- grid$Y
+ views <- list()
+
+ if ((k * tau) >= t_view[2]) {
+ N_Meanviews <- 1 #Number of views on expectations
+ v_mu_tmp <- array(0, dim = c(N_Meanviews, n_, T_))
+ v_mu_tmp[1, i_view[1], T_] <- 1
+ mu_view <- view[1]
+ } else{
+ N_Meanviews <- 2 #Number of views on expectations
+ v_mu_tmp <- array(0, dim = c(N_Meanviews, n_, T_))
+ v_mu_tmp[1, i_view[1], T_] <- 1
+ v_mu_tmp[2, i_view[2], T_ - 1] <- 1
+ mu_view <- view[1]
+ mu_view <- c( mu_view, view[2])
+ }
+
+ v_mu <- matrix(v_mu_tmp, nrow = nrow(v_mu_tmp))
+ views$N_Meanviews <- N_Meanviews
+ views$N_Covviews <- c()
+ views$dimension <- array(N)
+ views$monitoring_time <- array(T)
+ views$v_mu <- v_mu
+ views$v_sig <- NaN
+ views$mu_view <- mu_view
+ views$sig2_view <- c()
+
+ ############################################################################
+ #update the Posterior moments of the process
+ Posterior <- MVOU_Posterior(Prior, views)
+ alpha_t[1:n_, k] <- Posterior$mean_cost[(n_ + 1):(2 * n_)]
+ beta_t[1:n_, 1:n_, k] <- Posterior$mean_lin[(n_ + 1):(2 * n_), 1:n_] -
+ diag(1, n_)
+ sig2_1 <- Posterior$cov[(n_ + 1):(2 * n_), (n_ + 1):(2 * n_)]
+ sig2_1 <- omega %*% sig2_1 %*% t(omega)
+
+ ############################################################################
+ #update matrix q_t
+ psi_bb <- drop(psi_t_bb[1:k_, 1:k_, k])
+ q_t[1:k_, 1:k_, k] <- gamma * sig2_1 + eta * c2 + exp(-lambda) * psi_bb
+ }
+
+ b_prior <- matrix(NaN, t_, 1 )
+ b_post <- matrix(NaN, t_, 1 )
+
+ #Reconstructing the optimal exposure on the simulated path x
+ for (t in 1:t_) {
+ if (t == 1) {
+ b_legacy_prior <- b_legacy
+ b_legacy_post <- b_legacy
+ } else {
+ b_legacy_prior <- b_prior[t - 1]
+ b_legacy_post <- b_post[t - 1]
+ }
+ # prior exposure
+ b_prior[t] <- solve(q_prior) %*% (omega %*% alpha_prior + omega %*%
+ beta_prior %*% x[t,] + eta * c2 %*% b_legacy_prior +
+ exp(-lambda) %*% psi_bx_prior %*% (alpha_prior + (beta_prior +
+ diag(1, n_)) %*% x[t,]) + exp(-lambda) * psi_b_prior)
+
+ #posterior exposure
+ q <- drop(q_t[1:k_, 1:k_, t])
+ alpha <- drop(alpha_t[1:n_, t])
+ beta <- drop(beta_t[1:n_, 1:n_, t])
+ psi_bx <- drop(psi_t_bx[1:k_, 1:n_, t])
+ psi_b <- drop(psi_t_b[1:k_, t])
+ l <- omega %*% alpha + omega %*% beta %*% x[t,] + eta * c2 *
+ b_legacy_post + exp(-lambda) %*% psi_bx %*% (alpha + (beta +
+ diag(1, n_)) %*% x[t,]) + exp(-lambda) * psi_b
+ b_post[t] <- solve(q) %*% l
+ }
+ return(list(prior = b_prior, post = b_post))
+}
Modified: pkg/Meucci/R/FullyFlexibleBayesNets.R
===================================================================
--- pkg/Meucci/R/FullyFlexibleBayesNets.R 2015-07-31 09:27:25 UTC (rev 3894)
+++ pkg/Meucci/R/FullyFlexibleBayesNets.R 2015-08-02 17:21:04 UTC (rev 3895)
@@ -21,7 +21,7 @@
# initialize parameters
A <- matrix(, nrow = 0, ncol = nrow(X))
b <- g <- matrix(, nrow = 0, ncol = 1)
-
+ c = c()
# for each view...
for (k in 1:length(View)) {
I_mrg <- (X[, 1] < Inf)
@@ -64,8 +64,9 @@
A <- rbind(A, New_A) # constraint for the conditional expectations...
b <- rbind(b, New_b)
g <- rbind(g, -log(1 - View[[k]]$c))
+ c <- rbind(c,I_mrg)
}
- return(list(A = A, b = b, g = g))
+ return(list(A = A, b = b, g = g, c=c))
}
#' tweak a matrix
@@ -78,70 +79,29 @@
#' @author Ram Ahluwalia \email{ram@@wingedfootcapital.com}
Tweak <- function(A, b, g) {
- library(matlab)
- library(limSolve)
-
K <- nrow(A)
J <- ncol(A)
- g_ <- rbind(g, zeros(J, 1))
+ g_ <- rbind(g, matrix(0, J, 1))
- Aeq_ <- cbind(zeros(1, K), ones(1, J))
+ Aeq_ <- cbind(matrix(0, 1, K), matrix(1, 1, J))
beq_ <- 1
- lb_ <- rbind(zeros(K, 1), zeros(J, 1))
- ub_ <- rbind(Inf * ones(K, 1), ones(J, 1))
+ lb_ <- rbind(matrix(0, K, 1), matrix(0, J, 1))
+ ub_ <- rbind(Inf * matrix(1, K, 1), matrix(1, J, 1))
- A_ <- cbind(-eye(K), A)
+ A_ <- cbind(-diag(1, K), A)
b_ <- b
# add lower-bound and upper-bound constraints
- A_ <- rbind(A_, -eye(ncol(A_)))
- b_ <- rbind(b_, zeros(ncol(A_), 1))
+ A_ <- rbind(A_, -diag(1, ncol(A_)))
+ b_ <- rbind(b_, matrix(0, ncol(A_), 1))
- x0 <- rep(1/ncol(Aeq_), ncol(Aeq_))
-
- # db_ = linprog(g_, A_, b_, Aeq_,beq_, lb_, ub_) # MATLAB version
- # matrix containing coefficients of equality constraints Ex=F
- # optimResult = linp(E = Aeq_,
- # vector containing the right-hand side of equality constraints
- # F = beq_,
- # matrix containint coefficients of the inequality constraints GX >= H
- # G = -1*A_,
- # vector containing the right-hand side of the inequality constraints
- # H = -1*b_,
- # vector containing the coefficients of the cost function
- # Cost = -1*g_,
- # ispos = FALSE)
-
- costFunction <- function(x) {
- matrix(x, nrow = 1) %*% matrix(-1*g_, ncol = 1)
- }
- gradient <- function(x) {
- -1*g_
- }
- optimResult <- optim(par = x0,
- fn = costFunction, # CHECK
- gr = gradient,
- method = "L-BFGS-B",
- lower = lb_,
- upper = ub_,
- hessian = FALSE)
-
- # library(linprog)
- # optimResult2 = solveLP(E = Aeq_,
- # vector containing the right-hand side of equality constraints
- # F = beq_,
- # matrix containint coefficients of the inequality constraints GX >= H
- # G = -1*A_,
- # vector containing the right-hand side of the inequality constraints
- # H = -1*b_,
- # vector containing the coefficients of the cost function
- # Cost = -1*g_,
- # ispos = FALSE)
-
- db_ <- optimResult$X
-
+ x0 <- rep(1 / ncol(Aeq_), ncol(Aeq_))
+ dim(A_)
+ dim(b_)
+ dim(g_)
+ db_ <- solveLP(g_,A_,b_,Aeq_,beq_,lb_,ub_)
db <- db_[1:K]
return(db)
@@ -173,7 +133,6 @@
#' @export
ComputeMoments <- function(X, p) {
- library(matlab)
N <- ncol(X)
m <- t(X) %*% p
Sm <- t(X) %*% (X * repmat(p, 1, N)) # repmat : repeats/tiles a matrix
Added: pkg/Meucci/R/MVOUPosterior.R
===================================================================
--- pkg/Meucci/R/MVOUPosterior.R (rev 0)
+++ pkg/Meucci/R/MVOUPosterior.R 2015-08-02 17:21:04 UTC (rev 3895)
@@ -0,0 +1,106 @@
+#' @title Computes the posterior conditional expectation & covariance matrix of
+#'
+#' @description This function computes the posterior conditional expectation
+#' and covariance matrix at different monitoring dates t1,t2,...,t_ of the process
+#' X_{t1,t2...,t_} =(X_t1,
+#' X_t2,
+#' .
+#' .
+#' X_t_)
+#' X_t follows a MVOU process: dX_t <- (-theta*X_t+mu)dt + sig*dB_t
+#' The Views are: E{v_mu*X}=mu_view and Cov{v_sig*X}=sig2_view
+#'
+#' @param Mom list list of risk drivesrs
+#' @param Views list list of views
+#'
+#' @return Posterior list of posterior distribution information
+#'
+#' @details The list Mom consists of.
+#' \itemize{
+#' \item{Mom$monitoring_time }{ monitoring times =t_*n_ x 1]}
+#' \item{Mom$dimension }{ labels of the risk drivers [t_*n_ x 1]}
+#' \item{Mom$cov }{ prior covariance matrix of X_\{t1,t2...t_\} [t_*n_ x t_*n_]}
+#' \item{Mom$mean }{ prior vector of the means of X_\{t1,t2...t_\} [t_*n_ x 1]}
+#' \item{Mom$mean_cost}{[t_*n_ x 1]}
+#' \item{Mom$mean_lin}{[t_*n_ x 1]}
+#' \item{}{Mom$mean_cost and Mom$mean_lin are such that Mom$mean_cost +
+#' Mom$mean_lin*x0 = Mom$mean}}
+#'
+#' while the list of Views has the elements:
+#' \itemize{
+#' \item{Views$N_MeanViews }{ Number of views on the expectations [scalar]}
+#' \item{Views$N_CovViews }{ Number of views on the covariance matrix [scalar]}
+#' \item{Views$dimension }{ labels of the risk drivers [t_*n_ x 1]}
+#' \item{Views$monitoring_time }{ monitoring times [t_*n_ x 1]}
+#' \item{Views$v_mu }{ matrix that qualifies the views on expectations
+#' [N_MeanViews x t_*n_]}
+#' \item{Views$v_sig }{ matrix that qualifies the views on covariance
+#' [N_CovViews x t_*n_]}
+#' \item{Views$mu_view }{ extent of the views on expectation [N_MeanViews x 1]}
+#' \item{Views$sig2_view }{ extent of the views on the covariance
+#' [N_CovViews x 1]}}
+#'
+#' And the returned Posterior distribution list includes the elements:
+#' \itemize{
+#' \item{Posterior$monitoring_time }{ monitoring times [t_*n_ x 1]}
+#' \item{Posterior$dimension }{ labels of the risk drivers [t_*n_ x 1]}
+#' \item{Posterior$cov }{ posterior covariance matrice of X_{t1,t2...t_}
+#' [t_*n_ x t_*n_]}
+#' \item{Posterior$mean }{ posterior vector of the means of X_{t1,t2...t_}
+#' [t_*n_ x 1]}
+#' \item{Posterior$mean_cost [t_*n_ x 1]}
+#' \item{Posterior$mean_lin [t_*n_ x 1]}}
+#' Posterior$mean_cost and Posterior$mean_lin are such that
+#' Posterior$mean_cost + Posterior$mean_lin*x0 = Posterior$mean
+#'
+#' @references
+#' A. Meucci - "Dynamic Portfolio Management with Views at Multiple Horizons"
+#' \url{http://symmys.com/node/831}. See Meucci script for "MVOU_Prior.m"
+#'
+#' @author Xavier Valls \email{xavievallspla@@gmail.com}
+#' @export
+
+MVOU_Posterior <- function(Mom, Views) {
+ n <- unique(Mom$dimension)
+ t <- unique(Mom$monitoring_time)
+ n_ <- length(n)
+ t_ <- length(t)
+ grid <- meshgrid(t, 1:n_)
+ T <- grid$X
+ N <- grid$Y
+ Posterior <- list()
+ Posterior$monitoring_time <- array(T)
+ Posterior$dimension <- array(N)
+ Posterior$mean <- matrix(NaN, t_ * n_, 1)
+ Posterior$cov <- matrix(NaN, n_ * t_, n_ * t_)
+ Posterior$mean_cost <- matrix(NaN, t_ * n_, 1)
+ Posterior$mean_lin <- matrix(NaN, t_ * n_, n_)
+
+ S2 <- Mom$cov
+ mu <- Mom$mean
+ v_mu <- Views$v_mu
+ v_sig <- Views$v_sig
+ mu_view <- Views$mu_view
+ sig2_view <- Views$sig2_view
+
+ if (all(is.nan(Views$v_mu))) {
+ Posterior$mean <- mu
+ }else {
+ v_mu_dag <- (S2 %*% t(v_mu)) / (v_mu %*% S2 %*% t(v_mu))[1]
+ P_orth <- v_mu_dag %*% v_mu
+ P <- diag(1, dim(P_orth)[1], dim(P_orth)[2]) - P_orth
+ Posterior$mean <- P %*% mu + P_orth %*% v_mu_dag %*% mu_view
+ Posterior$mean_lin <- P %*% Mom$mean_lin
+ Posterior$mean_cost <- P %*% Mom$mean_cost + P_orth %*% v_mu_dag %*% mu_view
+ }
+ if (all(is.nan(Views$v_sig))) {
+ Posterior$cov <- S2
+ } else {
+ v_sig_dag <- (S2 %*% t(v_sig)) / (v_sig %*% S2 %*% t(v_sig))[1]
+ P_orth <- v_sig_dag %*% v_sig
+ diag(1, dim(P_orth)[1], dim(P_orth)[2]) - P_orth
+ Posterior$cov <- P %*% S2 %*% t(P) + P_orth %*% (v_sig_dag %*% sig2_view %*%
+ t(v_sig_dag)) %*% t(P_orth)
+ }
+ return(Posterior)
+}
Added: pkg/Meucci/R/MVOUPrior.R
===================================================================
--- pkg/Meucci/R/MVOUPrior.R (rev 0)
+++ pkg/Meucci/R/MVOUPrior.R 2015-08-02 17:21:04 UTC (rev 3895)
@@ -0,0 +1,106 @@
+#' @title Computes the conditional mean and covariance matrix at different
+#' dates
+#'
+#' @description This function computes the conditional mean and covariance
+#' matrix at different monitoring dates t1,t2,...,t_ of the process
+#' X_{t1,t2...,t_} =(X_t1,
+#' X_t2,
+#' .
+#' .
+#' X_t_)
+#' X_t follows a MVOU process: dX_t <- (-theta*X_t+mu)dt + sig*dB_t
+#'
+#' @param t [t_ x 1] vector of monitoring dates
+#' @param x0 [n_ x 1] observation at time 0
+#' @param theta [n_ x n_] transition matrix
+#' @param sig2 [n_ x n_] covariance matrix
+#' @param mu [n_ x 1] vector of drift parameters
+#'
+#' @return Mom list of the risk drivers
+#'
+#' @references
+#' A. Meucci - "Dynamic Portfolio Management with Views at Multiple Horizons"
+#' \url{http://symmys.com/node/831}. See Meucci script for "MVOU_Prior.m"
+#'
+#' @author Xavier Valls \email{xavievallspla@@gmail.com}
+#'
+#' @export
+
+#OUTPUT
+#Mom.monitoring_time <- monitoring times [t_*n_ x 1]
+#Mom.dimension <- labels of the risk drivers [t_*n_ x 1]
+#Mom.cov <- covariance matrix of X_{t1,t2...t_} [t_*n_ x t_*n_]
+#Mom.mean <- vector of the means of X_{t1,t2...t_} [t_*n_ x 1]
+#Mom.mean_cost [t_*n_ x 1]
+#Mom.mean_lin [t_*n_ x 1]
+#Mom.mean_cost and Mom.mean_lin are such that
+#t_ <- length(x)
+
+MVOU_Prior <- function (t, x0, theta, sig2, mu) {
+
+ Tol_eigb <- 10 ^ -8
+ t_ <- length(t)
+ n_ <- length(x0)
+ grid <- meshgrid(t, 1:n_)
+ Mom <- list()
+ Mom$monitoring_time <- grid$X
+ Mom$dimension <- t(grid$Y)
+ Mom$mean <- array(NaN, dim <- c(t_ * n_, 1))
+ Mom$cov <- array(NaN, dim <- c(t_ * n_, n_ * t_))
+ Mom$mean_cost <- array(NaN, dim <- c(t_ * n_, 1))
+ Mom$mean_lin <- array(NaN, dim <- c(t_ * n_, n_))
+
+ kronsum <- kronecker(theta, diag(1, n_)) + kronecker(diag(1, n_), theta)
+ e <- eigen(kronsum)
+ V <- e$vectors
+ D <- e$values
+ lambda_A <- array( NaN, length(D))
+
+ M <- matrix(NaN, n_, t_)
+ mean_lin <- array(NaN, dim <- c(n_, n_, t_))
+ mean_cost <- matrix(NaN, n_, t_)
+ e <- eigen(theta)
+ V1 <- e$vectors
+ theta_diag <- e$values
+
+ F <- array(NaN, n_)
+
+ for (i in 1:t_) {
+ F[theta_diag <= Tol_eigb] <- t[i]
+ F[theta_diag > Tol_eigb] <- (1 - exp(-theta_diag[theta_diag > Tol_eigb] *
+ t[i])) / theta_diag[theta_diag > Tol_eigb]
+ E <- expm(-theta * t[i])
+ M[ 1:n_, i] <- (E %*% x0 + V1 %*% diag(F) %*% pinv(V1) %*% mu)
+ M[, i] <- Re(M[,i])
+
+ mean_lin[1:n_, 1:n_, i] <- t(E)
+ mean_lin[ , , i] <- Re(mean_lin[, , i] )
+ mean_cost[1:n_, i] <- (V1 %*% diag(F) %*% pinv(V1) %*% mu)
+ mean_cost[ , i] <- Re(mean_cost[, i])
+
+ vecsig2 <- matrix(sig2, nrow = n_ ^ 2)
+ lambda_A[(abs(lambda) <= Tol_eigb)] <- t[i]
+ index <- abs(lambda) > Tol_eigb
+ lambda_A[index] <- (1 - exp(-lambda[index] * t[i])) / lambda[index]
+ A <- V %*% diag(lambda_A) %*% solve(V)
+ vecsig2_t <- A %*% vecsig2
+ sig2_t <- matrix(vecsig2_t, nrow = n_)
+ sig2_t <- Re(sig2_t)
+
+ for ( j in i:t_ ) {
+ Mom$cov[(i - 1) * n_ + (1:n_), (j - 1) * n_ + (1:n_)] <- sig2_t %*%
+ expm(-t(theta) *
+ (t[j] - t[i]))
+ Mom$cov[(j - 1) * n_ + (1:n_), (i - 1) * n_ + (1:n_)] <- expm(-theta *
+ (t[j] - t[i])) %*%
+ sig2_t
+ }
+ }
+ #Note that Mom.mean <- Mom.mean_cost + Mom.mean_lin*x0
+ Mom$mean <- matrix(M, ncol = 1)
+ Mom$mean_lin <- t(array(mean_lin, c(dim(mean_lin)[1], dim(mean_lin)[2] *
+ dim(mean_lin)[3])))
+ Mom$mean_cost <- t(matrix(mean_cost, nrow = 1))
+
+ return(Mom)
+}
Added: pkg/Meucci/R/QuadraticMatVC.R
===================================================================
--- pkg/Meucci/R/QuadraticMatVC.R (rev 0)
+++ pkg/Meucci/R/QuadraticMatVC.R 2015-08-02 17:21:04 UTC (rev 3895)
@@ -0,0 +1,75 @@
+#' @title Computes the matrix q_t of the problem to solve when using
+#' CALCULUS of VARIATION
+#'
+#' @description This function computes the matrix q_t of the problem to solve
+#' when using CALCULUS of VARIATION:
+#' argmin_b (b' q_t b - b'l_t)
+#'
+#' @param lambda [scalar] discounting parameter
+#' @param gamma [scalar] risk aversion parameter
+#' @param eta [scalar] overall weight of the market impact of transactions
+#' @param sig2 [n_*t_ x n_*t_] covariance matrix of the process of the risk
+#' drivers
+#' @param c2 [k_ x k_] market impact matrix
+#' @param tau_ [scalar] effective number of future time steps considered
+#' @param n_ [scalar] number of risk drivers
+#' @param i_invest [k_ x 1] labels of the investible risk drivers
+#'
+#' @return q_t [k_*tau_ x k_*tau_] matrix computed
+#' @note
+#' t_ = number of monitoring times at which sig2 is computed
+#' k_ = number of investible risk drivers
+#'
+#'
+#' @references
+#' A. Meucci - "Dynamic Portfolio Management with Views at Multiple Horizons"
+#' \url{http://symmys.com/node/831}. See Meucci script for "MVOU_Prior.m"
+#'
+#' @author Xavier Valls \email{xavievallspla@@gmail.com}
+#' @export
+
+QuadraticMat_Vc <- function(lambda, gamma, eta, sig2, c2, tau_, n_,
+ i_invest = NULL) {
+
+ if(length(i_invest) == 0)
+ i_invest <- 1:length(c2)
+
+
+ ExtractBlockMtx <- function( A, t1, t2, n1_, n2_, i_invest){
+ # This function extracts the (t1,t2)-block out of the block-diagonal matrix
+ # A and then it selects the indices given by i_invest
+ # A is a matrix t_*n1_ x t_*n2_. The matrix has t_ blocks. Each block is
+ # n1_ x n2_
+ blk <- A[(1 + (t1 - 1) * n1_):(t1 * n1_), (1 + (t2 - 1) * n2_):(t2 * n2_)]
+ if(length(blk) > 1)
+ blk <- blk[i_invest, i_invest]
+ return(blk)
+ }
+
+ k_ <- length(i_invest) #number of investible risk drivers
+ q_t <- matrix(0, k_ * tau_, k_ * tau_)
+ for (t in 1:(tau_ - 1)) {
+ sig2_t <- ExtractBlockMtx(sig2, t, t, n_, n_, i_invest)
+ sig2_t1 <- ExtractBlockMtx(sig2, t + 1, t + 1, n_, n_, i_invest)
+ sig2_tt1 <- ExtractBlockMtx(sig2, t, t + 1, n_, n_, i_invest)
+ sig2_t <- sig2_t + sig2_t1 - 2 * sig2_tt1
+ q_t[((t - 1) * k_ + 1):(t * k_), ((t - 1) * k_ + 1):(t * k_)] <-
+ exp(-lambda * (t - 1)) * (-gamma * 0.5 * sig2_t - eta * 0.5 * c2 *
+ (1 + exp(-lambda)))
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/returnanalytics -r 3895
More information about the Returnanalytics-commits
mailing list