[Returnanalytics-commits] r3745 - in pkg/Meucci: . R data demo man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu Jun 25 10:43:25 CEST 2015
Author: xavierv
Date: 2015-06-25 10:43:25 +0200 (Thu, 25 Jun 2015)
New Revision: 3745
Added:
pkg/Meucci/R/EffectiveBets.R
pkg/Meucci/R/InverseCallTransformation.R
pkg/Meucci/R/Torsion.R
pkg/Meucci/data/JGB.rda
pkg/Meucci/data/ghqx.rda
pkg/Meucci/data/linRet.rda
pkg/Meucci/data/pseudodata.rda
pkg/Meucci/demo/S_AnalyzeJGBrates.R
pkg/Meucci/demo/S_MainDiversification.R
pkg/Meucci/man/EffectiveBets.Rd
pkg/Meucci/man/InverseCallTransformation.Rd
pkg/Meucci/man/JGB.Rd
pkg/Meucci/man/PHist.Rd
pkg/Meucci/man/Torsion.Rd
pkg/Meucci/man/ghqx.Rd
pkg/Meucci/man/linRet.Rd
pkg/Meucci/man/pseudodata.Rd
Removed:
pkg/Meucci/man/pHist.Rd
Modified:
pkg/Meucci/NAMESPACE
pkg/Meucci/R/PlotDistributions.R
pkg/Meucci/R/RankingInformationFunctions.R
pkg/Meucci/demo/InvariantProjection.R
pkg/Meucci/demo/MeanDiversificationFrontier.R
pkg/Meucci/demo/Prior2Posterior.R
pkg/Meucci/demo/RankingInformation.R
pkg/Meucci/demo/logToArithmeticCovariance.R
pkg/Meucci/man/PlotFrontier.Rd
pkg/Meucci/man/PlotResults.Rd
pkg/Meucci/man/RIEfficientFrontier.Rd
pkg/Meucci/man/ViewRanking.Rd
Log:
- Added files missing. Fixed and formatted scripts up to RankingInformation demo
Modified: pkg/Meucci/NAMESPACE
===================================================================
--- pkg/Meucci/NAMESPACE 2015-06-24 22:54:29 UTC (rev 3744)
+++ pkg/Meucci/NAMESPACE 2015-06-25 08:43:25 UTC (rev 3745)
@@ -51,6 +51,7 @@
export(PlotDistributions)
export(PlotFrontier)
export(PlotMarginalsNormalInverseWishart)
+export(PlotResults)
export(PlotVolVsCompositionEfficientFrontier)
export(Prior2Posterior)
export(ProjectionStudentT)
@@ -68,6 +69,7 @@
export(TwoDimEllipsoid)
export(ViewCurveSlope)
export(ViewImpliedVol)
+export(ViewRanking)
export(ViewRealizedVol)
export(garch1f4)
export(garch2f8)
Added: pkg/Meucci/R/EffectiveBets.R
===================================================================
--- pkg/Meucci/R/EffectiveBets.R (rev 0)
+++ pkg/Meucci/R/EffectiveBets.R 2015-06-25 08:43:25 UTC (rev 3745)
@@ -0,0 +1,29 @@
+#' @title computes the Effective Number of Bets and the Diversification
+#' distribution
+#'
+#' @description computes the Effective Number of Bets and the Diversification
+#' distribution, as described in A. Meucci, A. Santangelo, R. Deguest -
+#' "Measuring Portfolio Diversification Based on Optimized Uncorrelated Factors"
+#'
+#' @param b [vector] (n_ x 1) exposures
+#' @param Sigma [matrix] (n_ x n_) covariance matrix
+#' @param t [matrix] (n_ x n_) torsion matrix
+#'
+#' @return enb [scalar] Effetive Number of Bets
+#' @return p [vector] (n_ x 1) diversification distribution
+#'
+#' @references
+#' A. Meucci, A. Santangelo, R. Deguest - "Measuring Portfolio Diversification
+#' Based on Optimized Uncorrelated Factors" \url{http://symmys.com/node/599}
+#'
+#' See Meucci's script "EffectiveBets.m"
+#
+#' @author Xavier Valls \email{xaviervallspla@@gmail.com}
+#' @export
+
+EffectiveBets <- function(b, Sigma, t) {
+ p <- mldivide(t(t), b) * (t %*% Sigma %*% b) / (t(b) %*% Sigma %*% b)[1];
+ enb <- exp(-sum(p * log(1 + (p - 1) * (p > 1e-5))));
+
+ return(list(p = p, enb = enb));
+}
Added: pkg/Meucci/R/InverseCallTransformation.R
===================================================================
--- pkg/Meucci/R/InverseCallTransformation.R (rev 0)
+++ pkg/Meucci/R/InverseCallTransformation.R 2015-06-25 08:43:25 UTC (rev 3745)
@@ -0,0 +1,61 @@
+#' @title Computes the Inverse Call Transformation
+#'
+#' @description Computes the Inverse Call Transformation and returns shadow
+#' rates, as described in A. Meucci, A. Loregian -"Neither Normal not Lognormal:
+#' Modeling Interest Rates Across all Regimes".
+#'
+#' @param Rates [matrix] (N x length(timeseries)) number of portfolio
+#' in the efficient frontier
+#' @param Tau [vector] (N x 1) vector containing the times to
+#' maturity corresponding to the rows of the rates matrix
+#' @param Eta [scalar] Inverse-call transformation parameter
+#' @param Zeta [scalar] Inverse-call transformation parameter
+#'
+#' @return x [vector] (NumPortf x 1) shadow rates, computed from
+#' rates via inverse-call transformation
+#'
+#' @note Smoothing parameter s obtained as s=eta*exp(zeta*tau)
+#'
+#' @references
+#' A. Meucci, A. Loregian - "Neither Normal not Lognormal: Modeling Interest
+#' Rates Across all Regimes" \url{http://symmys.com/node/601}
+#'
+#' See Meucci's script "InverseCallTransformation.m".
+#
+#' @author Xavier Valls \email{xaviervallspla@@gmail.com}
+#' @export
+
+InverseCallTransformation <- function(rates, tau, eta, zeta){
+
+ t_ <- dim(rates)[2]
+ x <- matrix(0, nrow(rates), ncol(rates))
+ s <- eta * exp( zeta * tau )
+
+
+ for (v in 1:length(tau)){
+ # inverse call transformation
+
+ x0 <- 0 # initialization
+
+ CallFit <- function(tmpX, y = rates[v, t], sigma = s[v]){
+ c <- BachelierCallPrice(tmpX, sigma)
+ return(F = y - c)
+ }
+ for (t in 1:t_){
+ # fitting inverse call
+ x[v,t] <- lsqnonlin(fun = CallFit, x0 = x0)$x
+ }
+ }
+
+ return(x)
+}
+
+# Bachelier call pricing function
+# Call function (zero-strike call option price profile according to the
+# Bachelier pricing function)
+#
+# s: smoothing parameter
+
+BachelierCallPrice <- function( x , s){
+ return(x * pnorm( x / s ) + s %*% dnorm(x / s))
+}
Modified: pkg/Meucci/R/PlotDistributions.R
===================================================================
--- pkg/Meucci/R/PlotDistributions.R 2015-06-24 22:54:29 UTC (rev 3744)
+++ pkg/Meucci/R/PlotDistributions.R 2015-06-25 08:43:25 UTC (rev 3745)
@@ -29,16 +29,16 @@
x <- as.matrix(seq(from = xl, to = xh, by = (xh - xl) / 100))
# posterior numerical
- PHist(X[,n] , p_ , NBins)
+ PHist(X[, n], p_ , NBins)
# posterior analytical
- y1 <- dnorm(x , Mu_[n] , sqrt(Sigma_[n,n]))
- lines(x , y1, type = "l", col = "red", xlab = "", ylab = "")
+ y1 <- dnorm(x, Mu_[n], sqrt(Sigma_[n,n]))
+ lines(x, y1, type = "l", col = "red", xlab = "", ylab = "")
# prior analytical
- y2 <- dnorm(x , Mu[n] ,sqrt(Sigma[n,n]))
- lines(x , y2, col = "blue", xlab = "", ylab = "")
+ y2 <- dnorm(x, Mu[n], sqrt(Sigma[n,n]))
+ lines(x, y2, col = "blue", xlab = "", ylab = "")
- legend(x = 1.5, y = 0.4 , legend = c("analytical","prior"),
- lwd = c(0.2,0.2), lty = c(1,1), col = c("red", "blue"))
+ legend(x = 1.5, y = 0.4 , legend = c("analytical", "prior"),
+ lwd = c(0.2, 0.2), lty = c(1, 1), col = c("red", "blue"))
}
-}
\ No newline at end of file
+}
Modified: pkg/Meucci/R/RankingInformationFunctions.R
===================================================================
--- pkg/Meucci/R/RankingInformationFunctions.R 2015-06-24 22:54:29 UTC (rev 3744)
+++ pkg/Meucci/R/RankingInformationFunctions.R 2015-06-25 08:43:25 UTC (rev 3745)
@@ -1,224 +1,271 @@
# TODO: add max weights constraint to EfficientFrontier()
# TODO: add computeCVaR to EfficientFrontier()
-# TODO: confirm QuadProg does not have a bug (i.e. it can optimize expected returns without use dvec by adding an equality constraint)
+# TODO: confirm QuadProg does not have a bug (i.e. it can optimize expected
+# returns without use dvec by adding an equality constraint)
-#' Plots the efficient frontier, as it appears in A. Meucci, "Fully Flexible Views: Theory and Practice", The Risk Magazine,
-#' October 2008, p 100-106.
+#' @title Plots the efficient frontier
#'
-#' @param e the NumPortf x 1 matrix of expected returns for each portfolio along the efficient frontier
-#' @param s the NumPortf x 1 matrix of standard deviation of returns for each portfolio along the efficient frontier
-#' @param w the NumPortf x N matrix of compositions (security weights) for each portfolio along the efficient frontier
-#'
-#' @references
-#' A. Meucci, "Fully Flexible Views: Theory and Practice" \url{http://www.symmys.com/node/158}
-#' See Meucci script for "RankingInformation/PlotFrontier.m"
+#' @description Plots the efficient frontier, as it appears in A. Meucci, "Fully
+#' Flexible Views: Theory and Practice", The Risk Magazine, October 2008,
+#' p 100-106.
#'
-#' @author Xavier Valls \email{flamejat@@gmail.com}
-#' @export
+#' @param e the NumPortf x 1 matrix of expected returns for each portfolio
+#' along the efficient frontier
+#' @param s the NumPortf x 1 matrix of standard deviation of returns for each
+#' portfolio along the efficient frontier
+#' @param w the NumPortf x N matrix of compositions (security weights) for
+#' each portfolio along the efficient frontier
+#'
+#' @references
+#' A. Meucci, "Fully Flexible Views: Theory and Practice"
+#' \url{http://www.symmys.com/node/158}. See Meucci script for
+#' "RankingInformation/PlotFrontier.m"
+#'
+#' @author Xavier Valls \email{xaviervallspla@@gmail.com}
+#' @export
-PlotFrontier = function( e, s, w )
-{
- xx = dim( w )[ 1 ];
- N = dim( w )[ 2 ];
- Data = t( apply( w, 1, cumsum ) );
+PlotFrontier <- function(e, s, w) {
+ N <- dim(w)[2]
+ Data <- t(apply(w, 1, cumsum))
- plot( c(min(s), 0), xlim = c( min(s) , max(s) ), ylim = c( 0, max(Data) ),
- main= "frontier", xlab = " Portfolio # risk propensity", ylab = "Portfolio composition" );
-
- for( n in 1 : N )
- {
- x = rbind( min(s), s, max(s) );
- y = rbind( 0, matrix( Data[ , N-n+1 ] ), 0 );
- polygon( x, y, col = rgb( 0.9 - mod(n,3)*0.2, 0.9 - mod(n,3)*0.2, 0.9 - mod(n,3)*0.2) );
+ plot(c(min(s), 0), xlim = c(min(s), max(s)), ylim = c(0, max(Data)),
+ main= "frontier", xlab = " Portfolio # risk propensity",
+ ylab = "Portfolio composition")
+
+ for (n in 1 : N) {
+ x <- rbind(min(s), s, max(s))
+ y <- rbind(0, matrix(Data[, N - n + 1]), 0)
+ polygon(x, y, col = rgb(0.9 - mod(n, 3) * 0.2, 0.9 - mod(n,3) * 0.2,
+ 0.9 - mod(n, 3) * 0.2))
}
}
-#' Plots the results of computing the efficient frontier (Expected returns and frontier), as it appears in A. Meucci, "Fully Flexible Views: Theory and Practice", The Risk Magazine,
-#' October 2008, p 100-106.
+#' @title Plots the results of computing the efficient frontier (Expected
+#' returns and frontier)
#'
-#' @param e the NumPortf x 1 matrix of expected returns for each portfolio along the efficient frontier
-#' @param s the NumPortf x 1 matrix of standard deviation of returns for each portfolio along the efficient frontier
-#' @param w the NumPortf x N matrix of compositions (security weights) for each portfolio along the efficient frontier
+#' @description Plots the results of computing the efficient frontier (Expected
+#' returns and frontier), as it appears in A. Meucci, "Fully Flexible Views:
+#' Theory and Practice", The Risk Magazine, October 2008, p 100-106.
+#'
+#' @param e the NumPortf x 1 matrix of expected returns for each portfolio
+#' along the efficient frontier
+#' @param s the NumPortf x 1 matrix of standard deviation of returns for
+#' each portfolio along the efficient frontier
+#' @param w the NumPortf x N matrix of compositions (security weights) for
+#' each portfolio along the efficient frontier
#' @param M the NumPortf x 1 vector of expected returns for each asset
#' @param Lower constraints
-#' @param Upper constraints
-#'
-#' @references
-#' A. Meucci, "Fully Flexible Views: Theory and Practice" \url{http://www.symmys.com/node/158}
-#' See Meucci script for "RankingInformation/PlotResults.m"
+#' @param Upper constraints
#'
-#' @author Xavier Valls \email{flamejat@@gmail.com}
+#' @references
+#' A. Meucci, "Fully Flexible Views: Theory and Practice"
+#' \url{http://www.symmys.com/node/158}. See Meucci script for
+#' "RankingInformation/PlotResults.m"
+#'
+#' @author Xavier Valls \email{xaviervallspla@@gmail.com}
+#' @export
-PlotResults = function( e, s, w, M, Lower = NULL , Upper = NULL)
-{
- N = length( M );
- dev.new();
- par( mfrow = c( 1, 2 ) );
- h1 = hist( M*100, plot = F )
- barplot( h1$density, horiz = T, main = "expected returns", xlab = "", ylab = "" );
- if(length(Lower) || length(Upper))
- {
- Changed = array( 0, N );
- Changed[ union( Lower, Upper ) ] = M[ union( Lower, Upper ) ] * 100;
- h2 = hist(Changed, plot = F );
- barplot( h2$density, horiz = T, col = "red", add = T );
+PlotResults <- function(e, s, w, M, Lower = NULL, Upper = NULL) {
+ N <- length(M)
+ dev.new()
+ par(mfrow = c(1, 2))
+ h1 <- hist(M * 100, plot = F)
+ barplot(h1$density, horiz = T, main = "expected returns", xlab = "",
+ ylab = "")
+ if (length(Lower) || length(Upper)) {
+ Changed <- array(0, N)
+ Changed[union(Lower, Upper)] <- M[union(Lower, Upper)] * 100
+ h2 <- hist(Changed, plot = F)
+ barplot(h2$density, horiz = T, col = "red", add = T)
}
- PlotFrontier( e*100, s*100, w );
+ PlotFrontier(e * 100, s * 100, w)
}
-
-#' Computes posterior probabilities to view the rankings, as it appears in A. Meucci,
-#' "Fully Flexible Views: Theory and Practice", The Risk Magazine, October 2008, p 100-106.
+#' @title Computes posterior probabilities to view the rankings
#'
+#' @description Computes posterior probabilities to view the rankings, as it
+#' appears in A. Meucci, "Fully Flexible Views: Theory and Practice", The Risk
+#' Magazine, October 2008, p 100-106.
+#'
#' @param X a vector containing returns for all the asset classes
#' @param p a vector containing the prior probability values
-#' @param Lower a vector of indexes indicating which column is lower than the corresponding column number in Upper
-#' @param Upper a vector of indexes indicating which column is lower than the corresponding column number in Upper
+#' @param Lower a vector of indexes indicating which column is lower than
+#' the corresponding column number in Upper
+#' @param Upper a vector of indexes indicating which column is lower than
+#' the corresponding column number in Upper
#'
#' @references
-#' A. Meucci, "Fully Flexible Views: Theory and Practice" \url{http://www.symmys.com/node/158}
-#' See Meucci script for "RankingInformation/ViewRanking.m"
+#' A. Meucci, "Fully Flexible Views: Theory and Practice"
+#' \url{http://www.symmys.com/node/158}. See Meucci script for
+#' "RankingInformation/ViewRanking.m"
#'
#' @author Ram Ahluwalia \email{ram@@wingedfootcapital.com}
-#' @export EntropyProg
+#' @export
-# example ViewRanking( X , p , Lower = c(3,4) , Upper = c(4,5) ) # two inequality views: asset 3 < asset 4 returns, and asset 4 < asset 5 returns
+# example ViewRanking(X, p, Lower = c(3,4), Upper = c(4,5))
+# two inequality views: asset 3 < asset 4 returns, and asset 4 < asset 5 returns
-ViewRanking = function( X , p , Lower , Upper )
-{
- library( matlab )
- J = nrow( X )
- N = ncol( X )
-
- K = length( Lower )
-
+ViewRanking <- function(X, p, Lower, Upper) {
+ J <- nrow(X)
+
# constrain probabilities to sum to one across all scenarios...
- Aeq = ones( 1 , J )
- beq = 1
-
+ Aeq <- matrix(1, 1, J)
+ beq <- 1
+
# ...constrain the expectations... A*x <= 0
# X[,Lower] refers to the column of returns for Asset-lower
# X[,Upper] refers to the column of returns for Asset-lower
- # X[ , Lower ] - X[ , Upper ] is vector returns of the "lower"" asset less the returns of the "higher" asset
- V = X[ , Lower ] - X[ , Upper ] # Jx1 vector. Expectation is assigned to each scenario
-
- A = t( V )
- b = 0 # The expectation is that (Lower - Upper)x <= 0. (i.e. The returns of upper are greater than zero for each scenario)
-
+ # X[, Lower] - X[, Upper] is vector returns of the "lower"" asset less the
+ # returns of the "higher" asset
+
+ # Jx1 vector. Expectation is assigned to each scenario
+ V <- X[, Lower] - X[, Upper]
+
+ A <- t(V)
+ # The expectation is that (Lower - Upper)x <= 0.
+ # (i.e. The returns of upper are greater than zero for each scenario)
+ b <- 0
+
# ...compute posterior probabilities
- p_ = EntropyProg( p , A , as.matrix(b) , Aeq , as.matrix(beq) )
-
- return( p_ )
+ p_ <- EntropyProg(p, A, as.matrix(b), Aeq, as.matrix(beq))
+
+ return(p_)
}
-#' Generates an efficient frontier based on Meucci's Ranking Information version and returns a A list with
-#' NumPortf efficient portfolios whos returns are equally spaced along the whole range of the efficient frontier,
-#' as it appears in A. Meucci, "Fully Flexible Views: Theory and Practice", The Risk Magazine, October 2008,
-#' p 100-106.
+#' @title Generates an efficient frontier
#'
+#' @description Generates an efficient frontier based on Meucci's Ranking
+#' Information version and returns a A list with NumPortf efficient portfolios
+#' whose returns are equally spaced along the whole range of the efficient
+#' frontier, as it appears in A. Meucci, "Fully Flexible Views: Theory and
+#' Practice", The Risk Magazine, October 2008, p 100-106.
+#'
#' Most recent version of article and MATLAB code available at
#' http://www.symmys.com/node/158
#'
-#' @param X a matrix with the joint-scenario probabilities by asset (rows are joint-scenarios, columns are assets)
-#' @param p a vector of probabilities associated with each scenario in matrix X
+#' @param X a matrix with the joint-scenario probabilities by asset
+#' (rows are joint-scenarios, columns are assets)
+#' @param p a vector of probabilities associated with each scenario
+#' in matrix X
#' @param Options a list of options....TBD
#'
-#' @return Exps the NumPortf x 1 vector of expected returns for each asset
-#' @return Covs the NumPortf x N vector of security volatilities along the efficient frontier
-#' @return w the NumPortf x N matrix of compositions (security weights) for each portfolio along the efficient frontier
-#' @return e the NumPortf x 1 matrix of expected returns for each portfolio along the efficient frontier
-#' @return s the NumPortf x 1 matrix of standard deviation of returns for each portfolio along the efficient frontier
+#' @return Exps the NumPortf x 1 vector of expected returns for each
+#' asset
+#' @return Covs the NumPortf x N vector of security volatilities along
+#' the efficient frontier
+#' @return w the NumPortf x N matrix of compositions (security
+#' weights)for each portfolio along the efficient frontier
+#' @return e the NumPortf x 1 matrix of expected returns for each
+#' portfolio along the efficient frontier
+#' @return s the NumPortf x 1 matrix of standard deviation of
+#' returns for each portfolio along the efficient frontier
#'
-#' @references
-#' A. Meucci, "Fully Flexible Views: Theory and Practice" \url{http://www.symmys.com/node/158}
-#' See Meucci script for "RankingInformation/EfficientFrontier.m"
+#' @references
+#' A. Meucci, "Fully Flexible Views: Theory and Practice"
+#' \url{http://www.symmys.com/node/158}. See Meucci script for
+#' "RankingInformation/EfficientFrontier.m"
#'
-#' @author Ram Ahluwalia \email{ram@@wingedfootcapital.com} and Xavier Valls \email{flamejat@@gmail.com}
+#' @author Ram Ahluwalia \email{ram@@wingedfootcapital.com} and
+#' Xavier Valls \email{xaviervallspla@@gmail.com}
#' @export
-RIEfficientFrontier = function( X , p , Options)
-{
+RIEfficientFrontier <- function(X, p, Options) {
- if( !require("limSolve") ) stop("This script requieres the limSolve package installed")
+ if(!require("limSolve"))
+ stop("This script requieres the limSolve package installed")
+ N <- ncol(X) # number of assets
- library( matlab )
-
- J = nrow( X ) # number of scenarios
- N = ncol( X ) # number of assets
-
- Exps = t(X) %*% p # probability-weighted expected return of each asset
-
- Scnd_Mom = t(X) %*% (X * ( p %*% matrix( 1, 1 , N ) ) )
- Scnd_Mom = ( Scnd_Mom + t(Scnd_Mom) ) / 2 # an N*N matrix
- Covs = Scnd_Mom - Exps %*% t( Exps )
-
- Constr = list()
-
+ Exps <- t(X) %*% p # probability-weighted expected return of each asset
+
+ Scnd_Mom <- t(X) %*% (X * (p %*% matrix(1, 1, N)))
+ Scnd_Mom <- (Scnd_Mom + t(Scnd_Mom)) / 2 # an N*N matrix
+ Covs <- Scnd_Mom - Exps %*% t(Exps)
+
+ Constr <- list()
+
# constrain the sum of weights to 1
- Constr$Aeq = matrix( 1, 1 , N )
- Constr$beq = 1
-
+ Constr$Aeq <- matrix(1, 1, N)
+ Constr$beq <- 1
+
# constrain the weight of any security to between 0 and 1
- Constr$Aleq = rbind( diag( 1, N ) , - diag( 1, N ) ) # linear coefficients matrix A in the inequality constraint A*x <= b
- Constr$bleq = rbind( matrix( 1, N, 1 ) , matrix( 0, N, 1 ) ) # constraint vector b in the inequality constraint A*x <= b
-
- Amat = rbind( Constr$Aeq , Constr$Aleq ) # stack the equality constraints on top of the inequality constraints
- bvec = rbind( Constr$beq , Constr$bleq ) # stack the equality constraints on top of the inequality constraints
-
- ############################################################################################
+
+ # linear coefficients matrix A in the inequality constraint A*x <= b
+ Constr$Aleq <- rbind(diag(1, N), - diag(1, N))
+ # constraint vector b in the inequality constraint A*x <= b
+ Constr$bleq <- rbind(matrix(1, N, 1), matrix(0, N, 1))
+
+ # stack the equality constraints on top of the inequality constraints
+ Amat <- rbind(Constr$Aeq, Constr$Aleq)
+ # stack the equality constraints on top of the inequality constraints
+ bvec <- rbind(Constr$beq, Constr$bleq)
+
+ ##############################################################################
# determine return of minimum-risk portfolio
- FirstDegree = matrix( 0, N , 1 ) # TODO: assumes that securities have zero expected returns when computing efficient frontier?
- SecondDegree = Covs
- # Why is FirstDegree "expected returns" set to 0?
+
+ # TODO: assumes that securities have zero expected returns when computing
+ # efficient frontier?
+ FirstDegree <- matrix(0, N, 1)
+ SecondDegree <- Covs
+ # Why is FirstDegree "expected returns" set to 0?
# We capture the equality view in the equality constraints matrix
- # In other words, we have a constraint that the Expected Returns by Asset %*% Weights = Target Return
- MinVol_Weights = solve.QP( Dmat = SecondDegree , dvec = -1*FirstDegree , Amat = -1*t(Amat) , bvec = -1*bvec , meq = length( Constr$beq ) )
- MinSDev_Exp = t( MinVol_Weights$solution ) %*% Exps
-
- ############################################################################################
+ # In other words, we have a constraint that the Expected Returns by
+ # Asset %*% Weights = Target Return
+ MinVol_Weights <- solve.QP(Dmat = SecondDegree, dvec = -1 * FirstDegree,
+ Amat = -1 * t(Amat), bvec = -1 * bvec,
+ meq = length(Constr$beq))
+ MinSDev_Exp <- t(MinVol_Weights$solution) %*% Exps
+
+ ##############################################################################
# determine return of maximum-return portfolio
- FirstDegree = -Exps
- MaxRet_Weights = linp( E = Constr$Aeq , F = Constr$beq , G = -1*Constr$Aleq , H = -1*Constr$bleq , Cost = FirstDegree , ispos = FALSE )$X
- MaxExp_Exp = t( MaxRet_Weights) %*% Exps
-
- ############################################################################################
+ FirstDegree <- -Exps
+ MaxRet_Weights <- linp(E = Constr$Aeq, F = Constr$beq, G = -1 * Constr$Aleq,
+ H = -1 * Constr$bleq, Cost = FirstDegree,
+ ispos = FALSE)$X
+ MaxExp_Exp <- t(MaxRet_Weights) %*% Exps
+
+ ##############################################################################
# slice efficient frontier in NumPortf equally thick horizontal sections
- Grid = matrix( , ncol = 0 , nrow = 0 )
- Grid = t( seq( from = Options$FrontierSpan[1] , to = Options$FrontierSpan[2] , length.out = Options$NumPortf ) )
-
- # the portfolio return varies from a minimum of MinSDev_Exp up to a maximum of MaxExp_Exp
- # We establish equally-spaced portfolio return targets and use this find efficient portfolios
- # in the next step
- Targets = as.numeric( MinSDev_Exp ) + Grid * as.numeric( ( MaxExp_Exp - MinSDev_Exp ) )
-
- ############################################################################################
- # compute the NumPortf compositions and risk-return coordinates
- FirstDegree = matrix( 0, N , 1 )
-
- w = matrix( , ncol = N , nrow = 0 )
- e = matrix( , ncol = 1 , nrow = 0 )
- s = matrix( , ncol = 1 , nrow = 0 )
-
- for ( i in 1:Options$NumPortf )
- {
+ Grid <- matrix(, ncol = 0, nrow = 0)
+ Grid <- t(seq(from = Options$FrontierSpan[1], to = Options$FrontierSpan[2],
+ length.out = Options$NumPortf))
+
+ # the portfolio return varies from a minimum of MinSDev_Exp up to a maximum of
+ # MaxExp_Exp. We establish equally-spaced portfolio return targets and use
+ # this find efficient portfolios in the next step.
+ Targets <- as.numeric(MinSDev_Exp) + Grid * as.numeric((MaxExp_Exp -
+ MinSDev_Exp))
+
+ ##############################################################################
+ # compute the NumPortf compositions and risk-return coordinates
+ FirstDegree <- matrix(0, N, 1)
+
+ w <- matrix(, ncol = N, nrow = 0)
+ e <- matrix(, ncol = 1, nrow = 0)
+ s <- matrix(, ncol = 1, nrow = 0)
+
+ for (i in 1:Options$NumPortf) {
# determine least risky portfolio for given expected return
- # Ax = b ; Exps %*% weights = Target Return
- AEq = rbind( Constr$Aeq , t( Exps ) ) # equality constraint: set expected return for each asset...
- bEq = rbind( Constr$beq , Targets[ i ] ) # ...and target portfolio return for i'th efficient portfolio
-
- Amat = rbind( AEq , Constr$Aleq ) # stack the equality constraints on top of the inequality constraints
- bvec = rbind( bEq , Constr$bleq )
-
- Weights = solve.QP( Dmat = SecondDegree , dvec = -1*FirstDegree , Amat = -1*t(Amat) , bvec = -1*bvec , meq = length( bEq ) )
-
- w = rbind( w , Weights$solution )
- s = rbind( s , sqrt( t(Weights$solution) %*% Covs %*% Weights$solution ) )
- e = rbind( e , Weights$solution %*% Exps )
+ # Ax = b Exps %*% weights = Target Return
+ # equality constraint: set expected return for each asset...
+ AEq <- rbind(Constr$Aeq, t(Exps))
+ # ...and target portfolio return for i'th efficient portfolio
+ bEq <- rbind(Constr$beq, Targets[i])
+
+ # stack the equality constraints on top of the inequality constraints
+ Amat <- rbind(AEq, Constr$Aleq)
+ bvec <- rbind(bEq, Constr$bleq)
+
+ Weights <- solve.QP(Dmat = SecondDegree, dvec = -1 * FirstDegree,
+ Amat = -1 * t(Amat), bvec = -1 * bvec, meq = length(bEq))
+
+ w <- rbind(w, Weights$solution)
+ s <- rbind(s, sqrt(t(Weights$solution) %*% Covs %*% Weights$solution))
+ e <- rbind(e, Weights$solution %*% Exps)
}
-
- return( list( e = e , Sdev = s , Composition = w , Exps = Exps , Covs = Covs ) )
+
+ return(list(e = e, Sdev = s, Composition = w, Exps = Exps, Covs = Covs))
}
Added: pkg/Meucci/R/Torsion.R
===================================================================
--- pkg/Meucci/R/Torsion.R (rev 0)
+++ pkg/Meucci/R/Torsion.R 2015-06-25 08:43:25 UTC (rev 3745)
@@ -0,0 +1,79 @@
+
+#' @title Computes the Principal Components torsion and the Minimum Torsion
+#'
+#' @description Computes the Principal Components torsion and the Minimum
+#' Torsion for diversification analysis, as described in A. Meucci,
+#' A. Santangelo, R. Deguest - "Measuring Portfolio Diversification Based on
+#' Optimized Uncorrelated Factors"
+#'
+#' @param Sigma [matrix] (n_ x n_) covariance matrix
+#' @param model [string] choose between 'pca' and 'minimum-torsion'
+#' model
+#' @param method [string] choose between 'approximate' and 'exact'
+#' method for 'minimum-torsion' model
+#' @param max_niter [scalar] choose number of iterations of numerical
+#' algorithm
+#'
+#' @return t [matrix] (n_ x n_) torsion matrix
+#'
+#' @references
+#' A. Meucci, A. Santangelo, R. Deguest - "Measuring Portfolio Diversification
+#' Based on Optimized Uncorrelated Factors" \url{http://symmys.com/node/599}
+#'
+#' See Meucci's script "torsion.m"
+#
+#' @author Xavier Valls \email{xaviervallspla@@gmail.com}
+#' @export
+
+Torsion = function(Sigma, model, method = NULL, max_niter = 10000 ) {
+
+ if (model == "pca") {
+ # PCA decomposition
+ e <- eigen(Sigma)
+ lambda <- e$values
+ ev <- e$vectors
+ flip <- ev[1,] < 0
+ # fix the sign of the eigenvector based on the sign of its first entry
+ ev[, flip] <- -ev[, flip]
+ index <- order(-lambda)
+
+ # PCA torsion
+ t <- t(ev[, index])
+
+ } else if (model == "minimum-torsion") {
+ # Correlation matrix
+ sigma <- diag(Sigma) ^ (1 / 2)
+ C <- diag(1 / sigma) %*% Sigma %*% diag(1 / sigma)
+ c <- sqrtm(C)$B # Riccati root of C
+
+ if (method == "approximate") {
+ t <- mrdivide(diag(sigma), c) %*% diag(1 / sigma)
+ } else if (method == "exact") {
+ n_ <- nrow(Sigma)
+
+ # initialize
+ d <- array(1, n_)
+ f <- array(0, max_niter)
+ for (i in 1:max_niter) {
+
+ U <- diag(d) %*% c %*% c %*% diag(d)
+ u <- sqrtm(U)$B
+ q <- mldivide(u, (diag(d) %*% c))
+ d <- diag(q %*% c)
+ pi_ <- diag(d) %*% q # perturbation
+ f[i] <- norm(c - pi_, type = "f")
+
+ if(i > 1 && abs(f[i] - f[i - 1]) / f[i] / n_ <= 10 ^ (-8)) {
+ f <- f[1:i]
+ break
+ } else if( i == max_niter && abs(f[max_niter] -
+ f[max_niter - 1]) / f[max_niter] / n_ > 10 ^ -8 ) {
+ print(paste("number of max iterations reached: n_iter =", max_niter))
+ }
+ }
+ x <- pi_ %*% solve(c)
+ t <- diag(sigma) %*% x %*% diag(1/sigma)
+ }
+ }
+ return(t)
+}
Added: pkg/Meucci/data/JGB.rda
===================================================================
(Binary files differ)
Property changes on: pkg/Meucci/data/JGB.rda
___________________________________________________________________
Added: svn:mime-type
+ application/octet-stream
Added: pkg/Meucci/data/ghqx.rda
===================================================================
(Binary files differ)
Property changes on: pkg/Meucci/data/ghqx.rda
___________________________________________________________________
Added: svn:mime-type
+ application/octet-stream
Added: pkg/Meucci/data/linRet.rda
===================================================================
(Binary files differ)
Property changes on: pkg/Meucci/data/linRet.rda
___________________________________________________________________
Added: svn:mime-type
+ application/octet-stream
Added: pkg/Meucci/data/pseudodata.rda
===================================================================
(Binary files differ)
Property changes on: pkg/Meucci/data/pseudodata.rda
___________________________________________________________________
Added: svn:mime-type
+ application/octet-stream
Modified: pkg/Meucci/demo/InvariantProjection.R
===================================================================
--- pkg/Meucci/demo/InvariantProjection.R 2015-06-24 22:54:29 UTC (rev 3744)
+++ pkg/Meucci/demo/InvariantProjection.R 2015-06-25 08:43:25 UTC (rev 3745)
@@ -1,65 +1,87 @@
# Annualization and Projection algorithm for invariant
#
-# SYMMYS - Last version of article and code available at http://symmys.com/node/136
+# SYMMYS - Last version of article and code available at
+# http://symmys.com/node/136
# Project summary statistics to arbitrary horizons under i.i.d. assumption
-# see Meucci, A. (2010) "Annualization and General Projection of Skewness, Kurtosis and All Summary Statistics"
+# see Meucci, A. (2010) "Annualization and General Projection of Skewness,
+# Kurtosis and All Summary Statistics"
# GARP Risk Professional, August, pp. 52-54
-N = 6 # a numeric with the number of the first N stadardized summary statistics to project
-K = 100 # a numeric with an arbitrary projection horizon
+# numeric with number of the first N stadardized summary statistics to project
+N <- 6
+# numeric with an arbitrary projection horizon
+K <- 100
-J = 100000 # a numeric with the number of scenarios
-a = -1 # a numeric with the location shift parameter. Mean of distribution will be exp(a)
-m = 0.2 # log of the mean of the distribution
-s = 0.4 # log of the standard deviation of the distribution
-
-X = GenerateLogNormalDistribution(J, a, m, s)
+# numeric with the number of scenarios
+J <- 100000
+# numeric with the location shift parameter. Mean of distribution will be exp(a)
+a <- -1
+# log of the mean of the distribution
+m <- 0.2
+# log of the standard deviation of the distribution
+s <- 0.4
+
+X <- GenerateLogNormalDistribution(J, a, m, s)
# TODO: Expectations on outputs
# Ga[1] should equal K*mean(X)
# Ga[2] should equal sqrt(K)*std(X)
-
-# show distribution of the invariant. Invariance test: The three distributions should be very similar
-hist( X , 50 , freq = FALSE , main = "Distribution of Invariant" , xlab = "X" ) # chart 1: distribution of invariant
-hist( X[ 1 : length( X ) / 2 ] , 50 , freq = FALSE , main = "Distribution (1st Half of Pop.)" , xlab = "X" ) # chart 2: distribution of invariant (1st-half of population)
-hist( X[ ( length( X ) / 2 ) : length( X ) ] , 50 , freq = FALSE , main = "Distribution (2nd Half of Pop.)" , xlab = "X" ) # chart 3: distribution of invariant (2nd-half of population)
-
+
+# show distribution of the invariant
+# Invariance test: The three distributions should be very similar
+# chart 1: distribution of invariant
+hist(X, 50, freq = FALSE, main = "Distribution of Invariant", xlab = "X")
+# chart 2: distribution of invariant (1st-half of population)
+hist(X[1 : length(X) / 2], 50, freq = FALSE,
+ main = "Distribution (1st Half of Pop.)", xlab = "X")
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/returnanalytics -r 3745
More information about the Returnanalytics-commits
mailing list