[Returnanalytics-commits] r2919 - in pkg/Meucci: . R demo man

noreply at r-forge.r-project.org noreply at r-forge.r-project.org
Wed Aug 28 20:12:16 CEST 2013


Author: xavierv
Date: 2013-08-28 20:12:15 +0200 (Wed, 28 Aug 2013)
New Revision: 2919

Added:
   pkg/Meucci/R/EfficientFrontierPrices.R
   pkg/Meucci/demo/S_MeanVarianceHorizon.R
   pkg/Meucci/demo/S_MeanVarianceOptimization.R
   pkg/Meucci/man/EfficientFrontierPrices.Rd
Modified:
   pkg/Meucci/DESCRIPTION
   pkg/Meucci/NAMESPACE
Log:
- added S_MeanVarianceHorizon and S_MeanVarianceOptimization demo scripts from chapter 6 and its associated functions

Modified: pkg/Meucci/DESCRIPTION
===================================================================
--- pkg/Meucci/DESCRIPTION	2013-08-28 18:02:06 UTC (rev 2918)
+++ pkg/Meucci/DESCRIPTION	2013-08-28 18:12:15 UTC (rev 2919)
@@ -92,5 +92,6 @@
     'MaxRsqCS.R'
     'EfficientFrontierReturnsBenchmark.R'
     'EfficientFrontierReturns.R'
+    'EfficientFrontierPrices.R'
     '
     FitOrnsteinUhlenbeck.R'

Modified: pkg/Meucci/NAMESPACE
===================================================================
--- pkg/Meucci/NAMESPACE	2013-08-28 18:02:06 UTC (rev 2918)
+++ pkg/Meucci/NAMESPACE	2013-08-28 18:12:15 UTC (rev 2919)
@@ -10,6 +10,7 @@
 export(ConvertCompoundedReturns2Price)
 export(Cumul2Raw)
 export(DetectOutliersViaMVE)
+export(EfficientFrontierPrices)
 export(EfficientFrontierReturns)
 export(EfficientFrontierReturnsBenchmark)
 export(EntropyProg)

Added: pkg/Meucci/R/EfficientFrontierPrices.R
===================================================================
--- pkg/Meucci/R/EfficientFrontierPrices.R	                        (rev 0)
+++ pkg/Meucci/R/EfficientFrontierPrices.R	2013-08-28 18:12:15 UTC (rev 2919)
@@ -0,0 +1,88 @@
+#' Compute the mean-variance efficient frontier (on prices) by quadratic programming, as described in 
+#' A. Meucci "Risk and Asset Allocation", Springer, 2005
+#'
+#'  @param   NumPortf       : [scalar] number of portfolio in the efficient frontier
+#'  @param   Covariance     : [matrix] (N x N) covariance matrix
+#'  @param   ExpectedValues : [vector] (N x 1) expected returns
+#'  @param   Current_Prices : [vector] (N x 1) current prices
+#'  @param   Budget         : [scalar] budget constraint
+#'  
+#'  @return  ExpectedValue  : [vector] (NumPortf x 1) expected values of the portfolios
+#'  @return  Std_Deviation  : [vector] (NumPortf x 1) standard deviations of the portfolios
+#'  @return  Composition    : [matrix] (NumPortf x N) optimal portfolios 
+#'
+#' @references
+#' \url{http://symmys.com/node/170}
+#' See Meucci's script for "EfficientFrontierReturns.m"
+#'
+#' @author Xavier Valls \email{flamejat@@gmail.com}
+#' @export
+
+EfficientFrontierPrices = function( NumPortf, Covariance, ExpectedValues, Current_Prices, Budget )
+{
+
+    NumAssets = ncol( Covariance );
+
+	##################################################################################################################
+	### Determine exp value of minimum-variance portfolio
+	FirstDegree  = matrix( 0, NumAssets, 1 );
+	SecondDegree = Covariance;
+	Aeq = t( Current_Prices );
+	beq = Budget;
+	A   = -diag( 1, NumAssets );
+	b   = matrix( 0, NumAssets, 1 );
+	Amat = rbind( Aeq, A);
+    bvec = rbind( beq, b);
+	#x0  = Budget / NumAssets * matrix( 1 , NumAssets, 1 );
+	MinVol_Allocation = matrix( solve.QP( Dmat = SecondDegree, dvec = -FirstDegree, Amat = -t(Amat), bvec = -bvec, meq = length( beq )  )$solution );
+	MinVol_ExpVal = t( MinVol_Allocation ) %*% ExpectedValues;
+
+	##################################################################################################################
+	### Determine exp value of maximum-expected value portfolio
+	Max_ExpVal = Budget * apply( ExpectedValues / Current_Prices, 2, max );
+
+	##################################################################################################################
+	### Slice efficient frontier in NumPortf equally thick horizontal sectors in the upper branch only
+	Target_ExpectedValues = MinVol_ExpVal + ( 0 : NumPortf ) * ( Max_ExpVal - MinVol_ExpVal ) / NumPortf;
+
+	##################################################################################################################
+	### Compute the NumPortf compositions and risk-return coordinates
+	Composition   = matrix( NaN, NumPortf, NumAssets );
+	Std_Deviation = matrix( NaN, NumPortf, 1 );
+	ExpectedValue = matrix( NaN, NumPortf, 1 );
+
+	Min_ExpectedValue = min(ExpectedValues);
+	Max_ExpectedValue = max(ExpectedValues);
+
+	IndexMin = which.min(ExpectedValues);	
+	IndexMax = which.max(ExpectedValues);
+	
+	for( i in 1 : NumPortf )
+	{
+	    # determine initial condition
+	    Matrix = rbind( cbind( Min_ExpectedValue, Max_ExpectedValue ),
+	              		cbind( Current_Prices[ IndexMin ], Current_Prices[ IndexMax ] ) );
+	          
+	    Allocation_0_MinMax = solve( Matrix ) %*% rbind( Target_ExpectedValues[ i ], Budget );
+	    
+	    Allocation_0 = matrix( 0, NumAssets, 1 );
+	    Allocation_0[ IndexMin ] = Allocation_0_MinMax[ 1 ];
+	    Allocation_0[ IndexMax ] = Allocation_0_MinMax[ 2 ];
+	    
+	    # determine least risky portfolio for given expected return
+	    AEq = rbind( Aeq, t(ExpectedValues) );
+	    bEq = rbind( beq, Target_ExpectedValues[ i ] );
+	    Amat = rbind( AEq, A);
+        bvec = rbind( bEq, b);
+
+	    #options = optimset('Algorithm', 'medium-scale');
+	    Allocation = t( solve.QP( Dmat = SecondDegree, dvec = -FirstDegree, Amat = -t(Amat), bvec = -bvec, meq = length( bEq )  )$solution );
+	    
+	    # store 
+	    Composition[ i, ]   = Allocation;
+	    Std_Deviation[ i ]  = sqrt(Allocation %*% Covariance %*% t(Allocation) );
+	    ExpectedValue[ i ]  = Target_ExpectedValues[ i ];
+	}
+
+	return( list( ExpectedValue = ExpectedValue, Std_Deviation = Std_Deviation, Composition = Composition ) );
+}
\ No newline at end of file

Added: pkg/Meucci/demo/S_MeanVarianceHorizon.R
===================================================================
--- pkg/Meucci/demo/S_MeanVarianceHorizon.R	                        (rev 0)
+++ pkg/Meucci/demo/S_MeanVarianceHorizon.R	2013-08-28 18:12:15 UTC (rev 2919)
@@ -0,0 +1,97 @@
+#' This script projects the distribution of the market invariants for stock market (i.e. compounded returns) 
+#' from the estimation interval to the investment horizon. 
+#' Then it computes the distribution of prices at the investment horizon and performs the two-step mean-variance
+#' optimization in terms of returns and relative portfolio weights.
+#' Described in A. Meucci,"Risk and Asset Allocation", Springer, 2005,  Chapter 6.
+#'
+#' @references
+#' \url{http://symmys.com/node/170}
+#' See Meucci's script for "S_MeanVarianceHorizon.m"
+#
+#' @author Xavier Valls \email{flamejat@@gmail.com}
+
+##################################################################################################################
+### Load data
+load("../data/stockSeries.Rda"); 	
+
+##################################################################################################################
+### Inputs
+tau = 1/252;        # time to horizon expressed in years
+tau_tilde = 1/52;  # estimation period expressed in years
+nSim   = 10000;
+Budget = 100;
+Zeta   = 10; # risk aversion parameter
+
+##################################################################################################################
+### Estimation of weekly invariants stock market (compounded returns)
+Week_C = diff( log( StockSeries$Prices.TimeSeries ) );
+N = ncol( Week_C );
+
+la = 0.1;
+Shrk_Exp = matrix( 0, N, 1);
+Exp_C_Hat = (1 - la) * matrix( apply( Week_C, 2, mean ) )  + la * Shrk_Exp;
+
+lb = 0.1;
+Shrk_Cov = diag( 1, N ) * sum( diag( cov( Week_C ) ) ) / N;
+Cov_C_Hat = (1 - lb) * cov(Week_C) + lb * (Shrk_Cov);
+
+##################################################################################################################
+### Stock market projection to horizon and pricing 
+Exp_Hrzn_C_Hat = Exp_C_Hat * tau / tau_tilde;
+Cov_Hrzn_C_Hat = Cov_C_Hat * tau / tau_tilde;
+StockCompReturns_Scenarios = rmvnorm( nSim, Exp_Hrzn_C_Hat, Cov_Hrzn_C_Hat);
+
+StockCurrent_Prices = matrix( StockSeries$Prices.TimeSeries[ nrow( StockSeries$Prices.TimeSeries ), ]);
+StockMarket_Scenarios = ( matrix( 1, nSim, 1) %*% t(StockCurrent_Prices)) * exp( StockCompReturns_Scenarios );
+##################################################################################################################
+### MV inputs - analytical
+Stock = ConvertCompoundedReturns2Price(Exp_Hrzn_C_Hat, Cov_Hrzn_C_Hat, StockCurrent_Prices);
+print( Stock$Exp_Prices );
+print( Stock$Cov_Prices );
+
+##################################################################################################################
+### MV inputs - numerical
+StockExp_Prices = matrix( apply( StockMarket_Scenarios, 2, mean ));
+StockCov_Prices = cov( StockMarket_Scenarios );
+print(StockExp_Prices);
+print(StockCov_Prices);
+
+StockExp_LinRets = StockExp_Prices / StockCurrent_Prices - 1;
+StockCov_LinRets = diag( c(1 / StockCurrent_Prices) ) %*% StockCov_Prices %*% diag( c(1 / StockCurrent_Prices) );
+
+##################################################################################################################
+### Portolio optimization
+# step 1: MV quadratic optimization to determine one-parameter frontier of quasi-optimal solutions ...
+NumPortf = 40;
+EFR = EfficientFrontierReturns( NumPortf, StockCov_LinRets, StockExp_LinRets );
+
+# step 2: ...evaluate satisfaction for all allocations on the frontier ...
+Store_Satisfaction = NULL;
+for( n in 1 : NumPortf )
+{
+  Allocation         = matrix( EFR$Composition[ n, ] ) * Budget / StockCurrent_Prices;
+  Objective_Scenario = StockMarket_Scenarios %*% Allocation;
+  Utility            = -exp( -1 / Zeta * Objective_Scenario);
+  ExpU               = apply( Utility, 2, mean );  
+  Satisfaction       = -Zeta * log( -ExpU );
+  Store_Satisfaction = cbind( Store_Satisfaction, Satisfaction ); ##ok<AGROW>
+}
+
+# ... and pick the best
+Optimal_Index 	   = which.max(Store_Satisfaction);
+Optimal_Allocation = EFR$Composition[ Optimal_Index, ];
+
+##################################################################################################################
+### Plots
+dev.new();
+
+par(mfrow = c( 2, 1 ) );
+# rets MV frontier 
+h = plot(EFR$Volatility, EFR$ExpectedValue, "l", lwd = 2, xlab = "st.dev. rets.", ylab = "exp.val rets.",
+ xlim = c( EFR$Volatility[1], EFR$Volatility[ length(EFR$Volatility) ]), ylim = c( min( EFR$ExpectedValue ), max( EFR$ExpectedValue ) ) );
+
+
+# satisfaction as function of st.deviation on the frontier
+h = plot( EFR$Volatility, Store_Satisfaction, "l", lwd = 2, xlab = "st.dev. rets.", ylab = "satisfaction",
+ xlim = c( EFR$Volatility[1], EFR$Volatility[ length(EFR$Volatility) ]), ylim = c( min(Store_Satisfaction), max(Store_Satisfaction) ) );
+

Added: pkg/Meucci/demo/S_MeanVarianceOptimization.R
===================================================================
--- pkg/Meucci/demo/S_MeanVarianceOptimization.R	                        (rev 0)
+++ pkg/Meucci/demo/S_MeanVarianceOptimization.R	2013-08-28 18:12:15 UTC (rev 2919)
@@ -0,0 +1,158 @@
+#' This script projects the distribution of the market invariants for the bond and stock markets 
+#' (i.e. the changes in yield to maturity and compounded returns) from the estimation interval to the investment 
+#' horizon
+#' Then it computes the distribution of prices at the investment horizon and performs the two-step mean-variance 
+#' optimization. 
+#' Described in A. Meucci,"Risk and Asset Allocation", Springer, 2005,  Chapter 6.
+#'
+#' @references
+#' \url{http://symmys.com/node/170}
+#' See Meucci's script for "S_MeanVarianceHorizon.m"
+#
+#' @author Xavier Valls \email{flamejat@@gmail.com}
+
+##################################################################################################################
+### Load data
+load( "../data/stockSeries.Rda" );
+
+##################################################################################################################
+### Inputs
+tau = 4 / 52;        # time to horizon expressed in years
+tau_tilde = 1 / 52;  # estimation period expressed in years
+
+FlatCurve  = 0.04;   
+TimesToMat = c( 4, 5, 10, 52, 520 ) / 52; # time to maturity of selected bonds expressed in years
+
+# parameters of the distribution of the changes in yield to maturity
+u_minus_tau = TimesToMat - tau;
+nu  = 8;
+mus = 0 * u_minus_tau;
+sigmas = ( 20 + 5 / 4 * u_minus_tau ) / 10000;
+
+Num_Scenarios = 100000;
+Budget = 100;
+Zeta = 10; # risk aversion parameter
+
+##################################################################################################################
+### Estimation of weekly invariants stock market (compounded returns)
+Week_C = diff( log( StockSeries$Prices.TimeSeries  ) );
+T = dim( Week_C )[1];
+N = dim( Week_C )[2]
+
+# shrinkage estimator of mean vector
+la = 0.1;
+Shrk_Exp = matrix( 0, N, 1 );
+Exp_C_Hat = (1 - la) * matrix( apply( Week_C, 2, mean) ) + la * Shrk_Exp;
+
+# shrinkage estimator of covariance
+lb = 0.1;
+Shrk_Cov = diag( 1, N ) * sum( diag( cov( Week_C ) ) ) / N;
+Cov_C_Hat = (1 - lb) * cov( Week_C ) + lb * ( Shrk_Cov );
+
+##################################################################################################################
+### Stock market projection to horizon and pricing 
+Exp_Hrzn_C_Hat = Exp_C_Hat * tau / tau_tilde;
+Cov_Hrzn_C_Hat = Cov_C_Hat * tau / tau_tilde;
+StockCompReturns_Scenarios = rmvnorm( Num_Scenarios, Exp_Hrzn_C_Hat, Cov_Hrzn_C_Hat );
+
+StockCurrent_Prices   = matrix( StockSeries$Prices.TimeSeries[ ncol(StockSeries$Prices.TimeSeries),  ] );
+StockMarket_Scenarios = ( matrix( 1, Num_Scenarios, 1 ) %*% t( StockCurrent_Prices ) ) * exp( StockCompReturns_Scenarios );
+
+##################################################################################################################
+### MV inputs
+# analytical
+CCR2P = ConvertCompoundedReturns2Price( Exp_Hrzn_C_Hat, Cov_Hrzn_C_Hat, StockCurrent_Prices ); 
+print( CCR2P$Exp_Prices );
+print( CCR2P$Cov_Prices );
+
+# numerical
+StockExp_Prices = matrix( apply( StockMarket_Scenarios, 2, mean) );
+StockCov_Prices = cov( StockMarket_Scenarios );
+print( StockExp_Prices );
+print( StockCov_Prices );
+
+##################################################################################################################
+### Bond market projection to horizon and pricing 
+BondCurrent_Prices_Shifted = exp( -FlatCurve * u_minus_tau );
+BondCurrent_Prices = exp( -FlatCurve * TimesToMat );
+
+# project bond market to horizon
+N = length( TimesToMat ); # number of bonds
+
+# generate common source of randomness
+U = runif( Num_Scenarios);
+BondMarket_Scenarios = matrix( 0, Num_Scenarios, N );
+for( n in 1 : N )
+{
+    # generate co-dependent changes in yield-to-maturity
+    DY_Scenarios = qnorm( U, mus[ n ] * tau / tau_tilde, sigmas[ n ] * sqrt( tau / tau_tilde ) ); 
+
+    # compute the horizon prices, (3.81) in "Risk and Asset Allocation" - Springer
+    X = -u_minus_tau[ n ] * DY_Scenarios;
+    BondMarket_Scenarios[ , n ] = BondCurrent_Prices_Shifted[ n ] * exp( X ); 
+}
+
+##################################################################################################################
+### MV inputs 
+
+# analytical
+Exp_Hrzn_DY_Hat  = mus * tau / tau_tilde;
+SDev_Hrzn_DY_Hat = sigmas * sqrt(tau / tau_tilde);
+Corr_Hrzn_DY_Hat = matrix( 1, N, N ); # full co-dependence
+Cov_Hrzn_DY_Hat  = diag(SDev_Hrzn_DY_Hat, length( SDev_Hrzn_DY_Hat)) %*% Corr_Hrzn_DY_Hat %*% diag(SDev_Hrzn_DY_Hat, length( SDev_Hrzn_DY_Hat));
+#[BondExp_Prices, BondCov_Prices]
+CCY2P = ConvertChangeInYield2Price(Exp_Hrzn_DY_Hat, Cov_Hrzn_DY_Hat, u_minus_tau, BondCurrent_Prices_Shifted);
+print( CCY2P$Exp_Prices );
+print( CCY2P$Cov_Prices );
+
+# numerical
+BondExp_Prices = matrix( apply( BondMarket_Scenarios,2, mean ) );
+BondCov_Prices = cov(BondMarket_Scenarios);
+print(BondExp_Prices);
+print(BondCov_Prices);
+
+##################################################################################################################
+### Portolio optimization
+# step 1: MV quadratic optimization to determine one-parameter frontier of quasi-optimal solutions ...
+E = rbind( StockExp_Prices, BondExp_Prices[ 2 ] );
+S = blkdiag( StockCov_Prices, matrix( BondCov_Prices[ 2, 2] ) );
+Current_Prices = rbind( StockCurrent_Prices, BondCurrent_Prices[ 2 ] );
+Market_Scenarios = cbind( StockMarket_Scenarios, BondMarket_Scenarios[ , 2 ] );
+
+NumPortf = 40;
+# frontier with QP (no short-sales)
+#[ExpectedValue, EFP$Std_Deviation, EFP$Composition] 
+
+EFP = EfficientFrontierPrices( NumPortf, S, E,Current_Prices, Budget );
+
+# step 2: ...evaluate satisfaction for all EFP$Composition on the frontier ...
+Store_Satisfaction = NULL;
+for( n in 1 : NumPortf )
+{
+  Allocation         = matrix( EFP$Composition[ n, ] );
+  Objective_Scenario = Market_Scenarios %*% Allocation;
+  Utility            = -exp( -1 / Zeta * Objective_Scenario );
+  ExpU               = apply( Utility, 2, mean );  
+  Satisfaction       = -Zeta * log( -ExpU );
+  Store_Satisfaction = cbind( Store_Satisfaction, Satisfaction ); ##ok<AGROW>
+}
+
+# ... and pick the best
+Optimal_Index = which.max( Store_Satisfaction );
+Optimal_Allocation = EFP$Composition[ Optimal_Index, ];
+print(Optimal_Allocation);
+
+##################################################################################################################
+### Plots
+dev.new()
+
+par(mfrow = c( 2, 1 ) );
+# rets MV frontier 
+h = plot( EFP$Std_Deviation, EFP$ExpectedValue, "l", lwd = 2, xlab = "st.dev. prices", ylab = "exp.val prices",
+ xlim = c( EFP$Std_Deviation[1], EFP$Std_Deviation[ length(EFP$Std_Deviation) ]), ylim = c( min(EFP$ExpectedValue), max(EFP$ExpectedValue) ) );
+
+
+# satisfaction as function of st.deviation on the frontier
+h = plot(EFP$Std_Deviation, Store_Satisfaction, "l", lwd = 2, xlab = "st.dev. prices", ylab = "satisfaction",
+ xlim = c( EFP$Std_Deviation[1], EFP$Std_Deviation[ length(EFP$Std_Deviation) ]), ylim = c( min(Store_Satisfaction), max(Store_Satisfaction) ) );
+

Added: pkg/Meucci/man/EfficientFrontierPrices.Rd
===================================================================
--- pkg/Meucci/man/EfficientFrontierPrices.Rd	                        (rev 0)
+++ pkg/Meucci/man/EfficientFrontierPrices.Rd	2013-08-28 18:12:15 UTC (rev 2919)
@@ -0,0 +1,43 @@
+\name{EfficientFrontierPrices}
+\alias{EfficientFrontierPrices}
+\title{Compute the mean-variance efficient frontier (on prices) by quadratic programming, as described in
+A. Meucci "Risk and Asset Allocation", Springer, 2005}
+\usage{
+  EfficientFrontierPrices(NumPortf, Covariance,
+    ExpectedValues, Current_Prices, Budget)
+}
+\arguments{
+  \item{NumPortf}{: [scalar] number of portfolio in the
+  efficient frontier}
+
+  \item{Covariance}{: [matrix] (N x N) covariance matrix}
+
+  \item{ExpectedValues}{: [vector] (N x 1) expected
+  returns}
+
+  \item{Current_Prices}{: [vector] (N x 1) current prices}
+
+  \item{Budget}{: [scalar] budget constraint}
+}
+\value{
+  ExpectedValue : [vector] (NumPortf x 1) expected values
+  of the portfolios
+
+  Std_Deviation : [vector] (NumPortf x 1) standard
+  deviations of the portfolios
+
+  Composition : [matrix] (NumPortf x N) optimal portfolios
+}
+\description{
+  Compute the mean-variance efficient frontier (on prices)
+  by quadratic programming, as described in A. Meucci "Risk
+  and Asset Allocation", Springer, 2005
+}
+\author{
+  Xavier Valls \email{flamejat at gmail.com}
+}
+\references{
+  \url{http://symmys.com/node/170} See Meucci's script for
+  "EfficientFrontierReturns.m"
+}
+



More information about the Returnanalytics-commits mailing list