[Returnanalytics-commits] r2582 - in pkg/Meucci: . R demo man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Tue Jul 16 17:31:33 CEST 2013
Author: xavierv
Date: 2013-07-16 17:31:33 +0200 (Tue, 16 Jul 2013)
New Revision: 2582
Added:
pkg/Meucci/R/GenerateUniformDrawsOnUnitSphere.R
pkg/Meucci/demo/S_BuyNHold.R
pkg/Meucci/demo/S_CPPI.R
pkg/Meucci/demo/S_CornishFisher.R
pkg/Meucci/demo/S_ESContributionFactors.R
pkg/Meucci/demo/S_ESContributionsStudentT.R
pkg/Meucci/demo/S_InvestorsObjective.R
pkg/Meucci/demo/S_UtilityMax.R
pkg/Meucci/demo/S_VaRContributionsUniform.R
pkg/Meucci/man/GenerateUniformDrawsOnUnitSphere.Rd
Modified:
pkg/Meucci/DESCRIPTION
pkg/Meucci/NAMESPACE
Log:
-added most of the demo scripts from chapter 5 and three from chapter 6
Modified: pkg/Meucci/DESCRIPTION
===================================================================
--- pkg/Meucci/DESCRIPTION 2013-07-16 04:59:26 UTC (rev 2581)
+++ pkg/Meucci/DESCRIPTION 2013-07-16 15:31:33 UTC (rev 2582)
@@ -82,3 +82,4 @@
'Raw2Cumul.R'
'FitExpectationMaximization.R'
'QuantileMixture.R'
+ 'GenerateUniformDrawsOnUnitSphere.R'
Modified: pkg/Meucci/NAMESPACE
===================================================================
--- pkg/Meucci/NAMESPACE 2013-07-16 04:59:26 UTC (rev 2581)
+++ pkg/Meucci/NAMESPACE 2013-07-16 15:31:33 UTC (rev 2582)
@@ -12,6 +12,7 @@
export(EntropyProg)
export(FitExpectationMaximization)
export(GenerateLogNormalDistribution)
+export(GenerateUniformDrawsOnUnitSphere)
export(hermitePolynomial)
export(integrateSubIntervals)
export(InterExtrapolate)
Added: pkg/Meucci/R/GenerateUniformDrawsOnUnitSphere.R
===================================================================
--- pkg/Meucci/R/GenerateUniformDrawsOnUnitSphere.R (rev 0)
+++ pkg/Meucci/R/GenerateUniformDrawsOnUnitSphere.R 2013-07-16 15:31:33 UTC (rev 2582)
@@ -0,0 +1,39 @@
+#' Generate a uniform sample on the unit hypersphere, as described in A. Meucci,
+#' "Risk and Asset Allocation", Springer, 2005.
+#'
+#' @param J : [scalar] number of draws
+#' @param N : [scalar] dimension
+#'
+#' @return X : [matrix] (T x N) of draws
+#'
+#'@note
+#' \item{ Initial script by Xiaoyu Wang - Dec 2006}
+#' \item{ We decompose X=U*R, where U is a uniform distribution on unit sphere and
+# R is a distribution on (0,1) proportional to r^(Dims-1), i.e. the area of surface of radius r }
+#'
+#' @references
+#' \url{http://symmys.com/node/170}
+#' See Meucci's script for "GenerateUniformDrawsOnUnitSphere.m"
+#'
+#' @author Xavier Valls \email{flamejat@@gmail.com}
+#' @export
+
+GenerateUniformDrawsOnUnitSphere = function(J, N)
+{
+l = matrix( 1, 1, N );
+
+# 1. Generate U
+Z = matrix( runif(J*N), J, N );
+normZ = sqrt( apply( Z * Z, 1, sum ) );
+U = Z / ( normZ %*% l );
+
+# 2. Generate R
+# the pdf of R is proportional to r^(N-1) therefore the cdf of R is r^N
+# we use quantile function of R sample R from uniform simulations
+Y = runif(J);
+R = Y ^ ( 1/N );
+
+# 3. Generate X
+X = U * ( R %*% l );
+ return( X );
+}
\ No newline at end of file
Added: pkg/Meucci/demo/S_BuyNHold.R
===================================================================
--- pkg/Meucci/demo/S_BuyNHold.R (rev 0)
+++ pkg/Meucci/demo/S_BuyNHold.R 2013-07-16 15:31:33 UTC (rev 2582)
@@ -0,0 +1,99 @@
+#' This script illustrates the buy & hold dynamic strategy, as 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_BuyNHold.m"
+#
+#' @author Xavier Valls \email{flamejat@@gmail.com}
+
+##################################################################################################################
+### Input parameters
+Initial_Investment = 1000;
+Time_Horizon = 6 / 12; # in years
+Time_Step = 1 / 252; # in years
+
+m = 0.2; # yearly expected return on the underlying
+s = 0.40; # yearly expected percentage volatility on the stock index
+r = 0.04; # risk-free (money market) interest rate
+
+NumSimul = 30000;
+
+##################################################################################################################
+# proportion of underlying you want to hold in the beginning, e.g.: 50
+Prct = 50 ;
+
+##################################################################################################################
+#### Initialize values
+Underlying_Index = Initial_Investment; # value of the underlyting at starting time, normalzed to equal investment
+Start = Underlying_Index;
+Elapsed_Time = 0;
+Portfolio_Value = Initial_Investment;
+
+Underlying_in_Portfolio_Percent = Prct / 100;
+
+Underlyings_in_Portfolio = Portfolio_Value * Underlying_in_Portfolio_Percent;
+Cash_in_Portfolio = Portfolio_Value - Underlyings_in_Portfolio;
+
+##################################################################################################################
+### Initialize parameters for the plot (no theory in this)
+
+Portfolio_Series = Portfolio_Value;
+Market_Series = Underlying_Index;
+Percentage_Series = Underlying_in_Portfolio_Percent;
+
+# asset evolution and portfolio rebalancing
+while( Elapsed_Time < (Time_Horizon - 10^(-5)) ) # add this term to avoid errors
+{
+ # time elapses...
+ Elapsed_Time = Elapsed_Time + Time_Step;
+
+ # ...asset prices evolve and portfolio takes on new value...
+ Multiplicator = exp( (m - s ^ 2 / 2) * Time_Step + s * sqrt( Time_Step ) * rnorm(NumSimul));
+ Underlying_Index = Underlying_Index * Multiplicator;
+ Underlyings_in_Portfolio = Underlyings_in_Portfolio * Multiplicator;
+ Cash_in_Portfolio = Cash_in_Portfolio * exp(r * Time_Step);
+ Portfolio_Value = Underlyings_in_Portfolio + Cash_in_Portfolio;
+
+ # ...and we rebalance our portfolio
+ Underlying_in_Portfolio_Percent = Underlyings_in_Portfolio / Portfolio_Value;
+
+ # store one path for the movie (no theory in this)
+ Portfolio_Series = cbind( Portfolio_Series, Portfolio_Value[ 1 ] ); ##ok<*AGROW>
+ Market_Series = cbind( Market_Series, Underlying_Index[ 1 ] );
+ Percentage_Series = cbind( Percentage_Series, Underlying_in_Portfolio_Percent[ 1 ] );
+}
+
+##################################################################################################################
+### Play the movie for one path
+Time = seq( 0, Time_Horizon, Time_Step);
+y_max = max( cbind( Portfolio_Series, Market_Series) ) * 1.2;
+dev.new();
+par( mfrow = c(2,1))
+for( i in 1 : length(Time) )
+{
+ plot( Time[ 1:i ], Portfolio_Series[ 1:i ], type ="l", lwd = 2.5, col = "blue", ylab = "value",
+ xlim = c(0, Time_Horizon), ylim = c(0, y_max), main = "investment (blue) vs underlying (red) value");
+ lines( Time[ 1:i ], Market_Series[ 1:i ], lwd = 2, col = "red" );
+ #axis( 1, [0, Time_Horizon, 0, y_max]);
+
+ plot(Time[ 1:i ], Percentage_Series[ 1:i ], type = "h", col = "red", xlab = "time", ylab = "#",
+ xlim = c(0, Time_Horizon), ylim =c(0,1), main = "percentage of underlying in portfolio");
+}
+
+##################################################################################################################
+### Plots
+# plot the scatterplot
+dev.new();
+
+# marginals
+NumBins = round(10 * log(NumSimul));
+layout( matrix(c(1,2,2,2,1,2,2,2,1,2,2,2,0,3,3,3), 4, 4, byrow = TRUE));
+barplot( table( cut( Portfolio_Value, NumBins )), horiz=TRUE, yaxt="n")
+
+# joint scatter plot
+plot(Underlying_Index, Portfolio_Value, xlab = "underlying at horizon (~ buy & hold )", ylab = "investment at horizon" );
+so = sort( Underlying_Index );
+lines( so, so, col = "red" );
+
+barplot( table( cut( Underlying_Index, NumBins )), yaxt="n")
Added: pkg/Meucci/demo/S_CPPI.R
===================================================================
--- pkg/Meucci/demo/S_CPPI.R (rev 0)
+++ pkg/Meucci/demo/S_CPPI.R 2013-07-16 15:31:33 UTC (rev 2582)
@@ -0,0 +1,112 @@
+#' This script illustrates the CPPI (constant proportion portfolio insurance) dynamic strategy, as 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_CPPI.m"
+#
+#' @author Xavier Valls \email{flamejat@@gmail.com}
+
+##################################################################################################################
+### Input parameters
+Initial_Investment = 1000;
+Time_Horizon = 6 / 12; # in years
+Time_Step = 1 / 252; # in years
+
+m = 0.2; # yearly expected return on the underlying
+s = 0.40; # yearly expected percentage volatility on the stock index
+r = 0.04; # risk-free (money market) interest rate
+
+nSim = 30000;
+
+##################################################################################################################
+### Setup
+# floor today (will evolve at the risk-free rate), e.g.: 950
+Floor = 980;
+# leverage on the cushion between your money and the floor, e.g. 3
+Multiple_CPPI = 5;
+
+##################################################################################################################
+### Initialize values
+Underlying_Index = Initial_Investment; # value of the underlyting at starting time, normalzed to equal investment
+Start = Underlying_Index;
+Elapsed_Time = 0;
+Portfolio_Value = Initial_Investment;
+
+Cushion = max( 0, Portfolio_Value - Floor );
+Underlyings_in_Portfolio = min(Portfolio_Value, max( 0, Multiple_CPPI * Cushion ) );
+#Cash_in_Portfolio = Portfolio_Value - Underlyings_in_Portfolio;
+Underlying_in_Portfolio_Percent = Underlyings_in_Portfolio / Portfolio_Value;
+
+Underlyings_in_Portfolio = Portfolio_Value * Underlying_in_Portfolio_Percent;
+Cash_in_Portfolio = Portfolio_Value - Underlyings_in_Portfolio;
+
+##################################################################################################################
+### Initialize parameters for the plot (no theory in this)
+Portfolio_Series = Portfolio_Value;
+Market_Series = Underlying_Index;
+Percentage_Series = Underlying_in_Portfolio_Percent;
+
+##################################################################################################################
+### Asset evolution and portfolio rebalancing
+while( Elapsed_Time < (Time_Horizon - 10^(-5)) ) # add this term to avoid errors
+{
+ # time elapses...
+ Elapsed_Time = Elapsed_Time + Time_Step;
+
+ # ...asset prices evolve and portfolio takes on new value...
+ Multiplicator = exp( (m - s ^ 2 / 2) * Time_Step + s * sqrt( Time_Step ) * rnorm( NumSimul ));
+ Underlying_Index = Underlying_Index * Multiplicator;
+ Underlyings_in_Portfolio = Underlyings_in_Portfolio * Multiplicator;
+ Cash_in_Portfolio = Cash_in_Portfolio * exp(r * Time_Step);
+ Portfolio_Value = Underlyings_in_Portfolio + Cash_in_Portfolio;
+
+ # ...and we rebalance our portfolio
+ Floor = Floor * exp( r * Time_Step );
+ Cushion = pmax( 0, (Portfolio_Value - Floor) );
+ Underlyings_in_Portfolio = pmin(Portfolio_Value, pmax( 0, Multiple_CPPI * Cushion) );
+ Cash_in_Portfolio = Portfolio_Value - Underlyings_in_Portfolio;
+ Underlying_in_Portfolio_Percent = Underlyings_in_Portfolio / Portfolio_Value;
+
+ # store one path for the movie (no theory in this)
+ Portfolio_Series = cbind( Portfolio_Series, Portfolio_Value[ 1 ] ); ##ok<*AGROW>
+ Market_Series = cbind( Market_Series, Underlying_Index[ 1 ] );
+ Percentage_Series = cbind( Percentage_Series, Underlying_in_Portfolio_Percent[ 1 ] );
+}
+
+##################################################################################################################
+### Play the movie for one path
+
+Time = seq( 0, Time_Horizon, Time_Step);
+y_max = max( cbind( Portfolio_Series, Market_Series) ) * 1.2;
+dev.new();
+par( mfrow = c(2,1))
+for( i in 1 : length(Time) )
+{
+ plot( Time[ 1:i ], Portfolio_Series[ 1:i ], type ="l", lwd = 2.5, col = "blue", ylab = "value",
+ xlim = c(0, Time_Horizon), ylim = c(0, y_max), main = "investment (blue) vs underlying (red) value");
+ lines( Time[ 1:i ], Market_Series[ 1:i ], lwd = 2, col = "red" );
+ #axis( 1, [0, Time_Horizon, 0, y_max]);
+
+ plot(Time[ 1:i ], Percentage_Series[ 1:i ], type = "h", col = "red", xlab = "time", ylab = "#",
+ xlim = c(0, Time_Horizon), ylim =c(0,1), main = "percentage of underlying in portfolio");
+
+}
+
+
+##################################################################################################################
+### plot the scatterplot
+dev.new();
+
+# marginals
+NumBins = round(10 * log(NumSimul));
+layout( matrix(c(1,2,2,2,1,2,2,2,1,2,2,2,0,3,3,3), 4, 4, byrow = TRUE));
+barplot( table( cut( Portfolio_Value, NumBins )), horiz=TRUE, yaxt="n");
+
+# joint scatter plot
+plot(Underlying_Index, Portfolio_Value, xlab = "underlying at horizon (~ buy & hold )", ylab = "investment at horizon" );
+so = sort( Underlying_Index );
+lines( so, so, col = "red" );
+
+barplot( table( cut( Underlying_Index, NumBins )), yaxt="n");
+
Added: pkg/Meucci/demo/S_CornishFisher.R
===================================================================
--- pkg/Meucci/demo/S_CornishFisher.R (rev 0)
+++ pkg/Meucci/demo/S_CornishFisher.R 2013-07-16 15:31:33 UTC (rev 2582)
@@ -0,0 +1,39 @@
+#'This script compares the Cornish-Fisher estiamte of the VaR with the true analytical VaR under the lognormal
+#'assumptions as described in A. Meucci,"Risk and Asset Allocation", Springer, 2005, Chapter 5.
+#'
+#' @references
+#' \url{http://symmys.com/node/170}
+#' See Meucci's script for "S_CornishFisher.m"
+#
+#' @author Xavier Valls \email{flamejat@@gmail.com}
+
+###################################################################################################################
+### Inputs
+mu = 0.05;
+sig = 0.05; # NB: change here and see the impact of approximation
+
+###################################################################################################################
+### Process data
+E_X = exp( mu + sig ^ 2 / 2 );
+Sd_X = exp( mu + sig ^ 2 / 2 ) * sqrt( exp( sig ^ 2 ) - 1 );
+Sk_X = sqrt( exp( sig ^ 2 ) - 1 ) * ( exp( sig ^ 2 ) + 2 );
+
+c = seq(0.001, 0.999, 0.001 );
+z = qnorm( c );
+
+Q_CF = E_X + Sd_X * ( z + Sk_X / 6 * ( z ^ 2 - 1 ) );
+Q_true = qlnorm( c,mu,sig );
+
+x = Q_true;
+f = dlnorm( x, mu, sig );
+
+###################################################################################################################
+### Plots
+dev.new();
+plot( x, f, type= "l", main = "pdf" );
+
+dev.new();
+plot( c, Q_true, type = "l", col = "red", main = "quantile" );
+lines( c, Q_CF );
+legend( "topleft", 1.9, c( "true", "Cornish-Fisher approx" ), col = c( "black","red" ),
+ lty=1, bg = "gray90" );
Added: pkg/Meucci/demo/S_ESContributionFactors.R
===================================================================
--- pkg/Meucci/demo/S_ESContributionFactors.R (rev 0)
+++ pkg/Meucci/demo/S_ESContributionFactors.R 2013-07-16 15:31:33 UTC (rev 2582)
@@ -0,0 +1,92 @@
+library(MASS);
+library(Matrix);
+#' This script computes the expected shortfall and the contributions to ES from each factor in simulations, using
+#' the conditional expectation definition of the contributions as described in A. Meucci,"Risk and Asset Allocation",
+#' Springer, 2005, Chapter 5.
+#'
+#' @references
+#' \url{http://symmys.com/node/170}
+#' See Meucci's script for "S_ESContributionFactors.m"
+#
+#' @author Xavier Valls \email{flamejat@@gmail.com}
+
+###################################################################################################################
+### Inputs
+
+N = 30; # number of securities
+K = 10; # number of factors
+a = runif(N); # allocation
+c = 0.95; # ES confidence
+
+###################################################################################################################
+### Generate market simulations
+# factor loadings
+B = matrix( runif( N*K ), N, K );
+
+# student t parameters
+nu = 10;
+eps = 0.1;
+g = 0.1;
+
+A = matrix( runif( K^2 ), K, K ) - 0.5;
+C = A %*% t(A);
+sigma_f = cov2cor(C);
+sigma_u = diag( 1, N );
+for( n in 1:(N - 1) )
+{
+ sigma_u = sigma_u + exp( -g * n ) *( rbind(cbind(matrix(0, N-n, n), diag(array( 1, N - n))), matrix(0,n,N)) +
+ rbind(matrix( 0, n, N ),cbind( diag( array( 1, N - n )), matrix( 0, N-n, n)) ));
+}
+
+sigma = as.matrix(.bdiag(list( eps * sigma_f, eps^2 * sigma_u))) #block diagonal matrix
+corr = cov2cor( sigma );
+diag_sigma = sqrt( diag( sigma ) );
+# scenarios
+nSim = 10000;
+l = matrix( 1, nSim);
+X = rmvt( nSim/2, corr, nu );
+X = rbind( X, -X ); # symmetrize simulations
+X = X %*% diag( diag_sigma );
+X = exp( X );
+F = X[ , 1:K ];
+U = X[ , (K+1):ncol(X) ];
+U = U - l %*% apply( U, 2, mean );
+M = F %*% t( B ) + U;
+
+###################################################################################################################
+### Risk management
+# compute the objective
+Psi = M %*% a;
+
+# compute ES
+th = ceiling((1-c) * nSim); # threshold
+spc = matrix( 0, nSim, 1 );
+spc[ 1 : th ] = 1;
+spc = spc / sum( spc );
+
+Sort_Psi = sort( Psi );
+Index = order( Psi );
+ES_simul = t(Sort_Psi) %*% spc;
+
+# augment factor set to include residual
+F_ = cbind( F, U%*%a );
+# compute portfolio-level loadings
+b_ = cbind( t(a)%*%B, 1 );
+
+# sort factors according to order induced by objective's realizations
+
+Sort_F_ = F_[Index, ];
+DES_simul = matrix( NaN, 1, K+1 );
+for( k in 1 : (K+1) )
+{
+ # compute gradient as conditional expectation
+ DES_simul[ k ] = t(spc) %*% Sort_F_[ , k ];
+}
+# compute contributions
+ContrES_simul = b_ * DES_simul;
+
+###################################################################################################################
+### Plots
+dev.new();
+bar = barplot( ContrES_simul, main = "contributions to ES" );
+axis( 1 , at = bar, labels= 1:ncol(ContrES_simul));
Added: pkg/Meucci/demo/S_ESContributionsStudentT.R
===================================================================
--- pkg/Meucci/demo/S_ESContributionsStudentT.R (rev 0)
+++ pkg/Meucci/demo/S_ESContributionsStudentT.R 2013-07-16 15:31:33 UTC (rev 2582)
@@ -0,0 +1,89 @@
+library(MASS);
+library(Matrix);
+#' This script computes the expected shortfall and the contributions to ES from each security:
+#' - analytically, under the Student t assumption for the market
+#' - in simulations, using the conditional expectation definition of the contributions
+#' Described in A. Meucci,"Risk and Asset Allocation",Springer, 2005, Chapter 5.
+#'
+#' @references
+#' \url{http://symmys.com/node/170}
+#' See Meucci's script for "S_ESContributionsStudentT.m"
+#
+#' @author Xavier Valls \email{flamejat@@gmail.com}
+
+##################################################################################################################
+### Inputs
+# number of assets
+N = 40;
+
+# market parameters (student t distribution)
+Mu = runif(N);
+A = matrix( runif( N^2 ), N, N )-.5;
+Sigma = A %*% t(A);
+nu = 7;
+
+# allocation
+a = runif(N) - 0.5;
+
+# ES confidence
+c = 0.95;
+
+nSim = 10000;
+
+###################################################################################################################
+### Generate market scenarios
+l = matrix( 1, nSim, 1);
+diag_s = diag( sqrt( diag( Sigma ) ) );
+invdiag_s = diag( 1 / sqrt( diag( Sigma ) ) );
+C = invdiag_s * Sigma * invdiag_s;
+X = rmvt( nSim / 2, C, nu);
+X = rbind( X, -X ); # symmetrize simulations
+M = l %*% t( Mu ) + X %*% diag_s;
+
+###################################################################################################################
+### Simulations
+# compute the objective
+Psi = M %*% a;
+
+# compute cut-off spectrum (step function) for empirical ES estimation, see (5.218)
+th = ceiling((1-c) * nSim); # threshold
+spc = matrix( 0, nSim, 1 );
+spc[ 1 : th ] = 1;
+spc = spc / sum( spc );
+
+# compute ES
+Sort_Psi = sort( Psi );
+Index = order( Psi );
+ES_simul = t(Sort_Psi) %*% spc;
+
+# sort market according to order induced by objective's realizations
+Sort_M = M[ Index, ];
+DES_simul = matrix(NaN, 1, N);
+for( n in 1 : N )
+{
+ # compute gradient as conditional expectation
+ DES_simul[ n ] = t(spc) %*% Sort_M[ , n ];
+}
+
+# compute contributions
+ContrES_simul = a * t( DES_simul);
+
+###################################################################################################################
+### Analytical
+# this does NOT depend on the allocation...
+ES_standardized = 1 / ( 1 - c ) * integrate(qt, lower=10^(-8) ,upper=(1-c), df=7)$value;
+
+# ...the dependence on the allocation is analytical
+ES_an = t( Mu ) %*% a + ES_standardized * sqrt(t(a) %*% Sigma %*% a);
+DES_an = Mu + ES_standardized * Sigma %*% a / sqrt( t(a) %*% Sigma %*% a)[1];
+ContrES_an = a * DES_an;
+
+###################################################################################################################
+### Plots
+dev.new();
+par( mfrow = c( 2, 1 ) );
+bar = barplot(t(ContrES_an), main = "contributions to ES, analytical" );
+axis( 1 , at = bar, labels= 1:length(ContrES_an[,1]));
+
+bar = barplot(t(ContrES_simul), main = "contributions to ES, simulations" );
+axis( 1 , at = bar, labels= 1:length(ContrES_simul[,1]));
\ No newline at end of file
Added: pkg/Meucci/demo/S_InvestorsObjective.R
===================================================================
--- pkg/Meucci/demo/S_InvestorsObjective.R (rev 0)
+++ pkg/Meucci/demo/S_InvestorsObjective.R 2013-07-16 15:31:33 UTC (rev 2582)
@@ -0,0 +1,79 @@
+library(mvtnorm);
+#' This script familiarizes the users with the objectives of different investors in a highly
+#' non-normal bi-variate market of securities, as described in A. Meucci,"Risk and Asset
+#' Allocation",Springer, 2005, Chapter 5.
+#'
+#' @references
+#' \url{http://symmys.com/node/170}
+#' See Meucci's script for "S_InvestorsObjective.m"
+#
+#' @author Xavier Valls \email{flamejat@@gmail.com}
+
+##################################################################################################################
+### Parameters of first marginal
+nu_1 = 3;
+s2_1 = 3;
+
+# parameters of second marginal
+mu_2 = 0.1;
+s2_2 = 0.2;
+
+# correlation in normal copula
+r = 0.5;
+
+# number of simulations
+J = 10000;
+
+# portfolio allocation
+a = matrix(c( 1, 2 ));
+
+# benchmark allocation
+b = matrix( c( 2, 1 ));
+
+##################################################################################################################
+### Compute current prices
+p_1 = nu_1 * s2_1;
+p_2 = exp( mu_2 + 0.5 * s2_2 ^ 2 );
+p = matrix( c( p_1, p_2 ));
+
+##################################################################################################################
+### Generate samnple of prices at the investment horizon
+N = rmvnorm(J, cbind( 0, 0 ), rbind( c(1, r), c(r, 1)));
+N_1 = N[ , 1 ];
+N_2 = N[ , 2 ];
+
+U_1 = pnorm( N_1 );
+U_2 = pnorm( N_2 );
+
+aa = nu_1 / 2;
+bb = 2 * s2_1;
+P_1 = qgamma( U_1, aa, scale = bb);
+P_2 = qlnorm( U_2, mu_2, sqrt(s2_2));
+
+P = cbind( P_1, P_2 );
+
+# generate sample of final wealth
+W = P %*% a;
+
+# generate sample of PnL
+PnL = (P - matrix( 1, J, 1) %*% t( p )) %*% a;
+
+# generate sample of benchmark-relative wealth
+K = diag(1, 2) - p %*% t(b) / (t(b) %*% p)[1];
+WRel = P %*% t(K) %*% a;
+
+##################################################################################################################
+### Plots
+NumBins = round(10 * log(J));
+dev.new();
+plot(P_1, P_2, xlab = "P_1", ylab = "P_2" );
+
+dev.new()
+hist(W, NumBins, main = "final wealth");
+
+dev.new();
+hist(PnL, NumBins, main = "P&L");
+
+dev.new();
+hist(WRel, NumBins, main = "benchmark-relative wealth" );
+
Added: pkg/Meucci/demo/S_UtilityMax.R
===================================================================
--- pkg/Meucci/demo/S_UtilityMax.R (rev 0)
+++ pkg/Meucci/demo/S_UtilityMax.R 2013-07-16 15:31:33 UTC (rev 2582)
@@ -0,0 +1,101 @@
+#' This script illustrates the constant weight dynamic strategy that maximizes power utility, as 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_UtilityMax.m"
+#
+#' @author Xavier Valls \email{flamejat@@gmail.com}
+
+##################################################################################################################
+### Input parameters
+Initial_Investment = 1000;
+Time_Horizon = 6/12; # in years
+Time_Step = 1/252; # in years
+
+m = 0.2; # yearly expected return on the underlying
+s = 0.40; # yearly expected percentage volatility on the stock index
+r = 0.04; # risk-free (money market) interest rate
+
+nSim = 30000;
+
+# amount to be invested in the underlying? e.g.: 50
+Prct = 50 ;
+
+##################################################################################################################
+### Initialize values
+Underlying_Index = Initial_Investment; # value of the underlyting at starting time, normalzed to equal investment
+Start = Underlying_Index;
+Elapsed_Time = 0;
+Portfolio_Value = Initial_Investment;
+
+Underlying_in_Portfolio_Percent = Prct / 100;
+
+Underlyings_in_Portfolio = Portfolio_Value * Underlying_in_Portfolio_Percent;
+Cash_in_Portfolio = Portfolio_Value - Underlyings_in_Portfolio;
+
+##################################################################################################################
+### Initialize parameters for the plot (no theory in this)
+Portfolio_Series = Portfolio_Value;
+Market_Series = Underlying_Index;
+Percentage_Series = Underlying_in_Portfolio_Percent;
+
+# asset evolution and portfolio rebalancing
+while( Elapsed_Time < (Time_Horizon - 10^(-5)) ) # add this term to avoid errors
+{
+ # time elapses...
+ Elapsed_Time = Elapsed_Time + Time_Step;
+
+ # ...asset prices evolve and portfolio takes on new value...
+ Multiplicator = exp( (m - s ^ 2 / 2) * Time_Step + s * sqrt( Time_Step ) * rnorm( NumSimul ));
+ Underlying_Index = Underlying_Index * Multiplicator;
+ Underlyings_in_Portfolio = Underlyings_in_Portfolio * Multiplicator;
+ Cash_in_Portfolio = Cash_in_Portfolio * exp(r * Time_Step);
+ Portfolio_Value = Underlyings_in_Portfolio + Cash_in_Portfolio;
+
+ # ...and we rebalance our portfolio
+ #Underlying_in_Portfolio_Percent = Underlying_in_Portfolio_Percent;
+ Underlyings_in_Portfolio = Portfolio_Value * Underlying_in_Portfolio_Percent;
+ Cash_in_Portfolio = Portfolio_Value - Underlyings_in_Portfolio;
+
+ # store one path for the movie (no theory in this)
+ Portfolio_Series = cbind( Portfolio_Series, Portfolio_Value[ 1 ] ); ##ok<*AGROW>
+ Market_Series = cbind( Market_Series, Underlying_Index[ 1 ] );
+ Percentage_Series = cbind( Percentage_Series, Underlying_in_Portfolio_Percent[ 1 ] );
+}
+
+##################################################################################################################
+### Play the movie for one path
+
+Time = seq( 0, Time_Horizon, Time_Step);
+y_max = max( cbind( Portfolio_Series, Market_Series) ) * 1.2;
+dev.new();
+par( mfrow = c(2,1))
+for( i in 1 : length(Time) )
+{
+ plot( Time[ 1:i ], Portfolio_Series[ 1:i ], type ="l", lwd = 2.5, col = "blue", ylab = "value",
+ xlim = c(0, Time_Horizon), ylim = c(0, y_max), main = "investment (blue) vs underlying (red) value");
+ lines( Time[ 1:i ], Market_Series[ 1:i ], lwd = 2, col = "red" );
+ #axis( 1, [0, Time_Horizon, 0, y_max]);
+
+ plot(Time[ 1:i ], Percentage_Series[ 1:i ], type = "h", col = "red", xlab = "time", ylab = "#",
+ xlim = c(0, Time_Horizon), ylim =c(0,1), main = "percentage of underlying in portfolio");
+
+}
+
+##################################################################################################################
+### Plots
+# plot the scatterplot
+dev.new();
+
+# marginals
+NumBins = round(10 * log(NumSimul));
+layout( matrix(c(1,2,2,2,1,2,2,2,1,2,2,2,0,3,3,3), 4, 4, byrow = TRUE));
+barplot( table( cut( Portfolio_Value, NumBins )), horiz=TRUE, yaxt="n");
+
+# joint scatter plot
+plot(Underlying_Index, Portfolio_Value, xlab = "underlying at horizon (~ buy & hold )", ylab = "investment at horizon" );
+so = sort( Underlying_Index );
+lines( so, so, col = "red" );
+
+barplot( table( cut( Underlying_Index, NumBins )), yaxt="n");
\ No newline at end of file
Added: pkg/Meucci/demo/S_VaRContributionsUniform.R
===================================================================
--- pkg/Meucci/demo/S_VaRContributionsUniform.R (rev 0)
+++ pkg/Meucci/demo/S_VaRContributionsUniform.R 2013-07-16 15:31:33 UTC (rev 2582)
@@ -0,0 +1,76 @@
+
+#' This script computes the VaR and the contributions to VaR from each security
+#' - analytically, under the elliptical-uniform assumption for the market
+#' - in simulations, using the conditional expectation definition of the contributions
+#' Described in A. Meucci,"Risk and Asset
+#' Allocation",Springer, 2005, Chapter 5.
+#'
+#' @references
+#' \url{http://symmys.com/node/170}
+#' See Meucci's script for "S_VaRContributionsUniform.m"
+#
+#' @author Xavier Valls \email{flamejat@@gmail.com}
+
+###################################################################################################################
+### Inputs
+# number of assets
+N = 10;
+
+# market parameters (uniform on ellipsoid)
+Mu = matrix(runif(N));
+A = matrix( runif(N*N), N, N) - 0.5;
+Sigma = A * t(A);
+
+# allocation
+a = matrix( runif(N) ) - 0.5;
+
+# quantile level
+c = 0.95;
+
+nSim = 10000;
+
+###################################################################################################################
+### Generate market scenarios
+X = GenerateUniformDrawsOnUnitSphere(nSim, N); # uniform on sphere
+M = matrix( 1, nSim, 1) %*% t(Mu) + X %*% t(A);
+
+###################################################################################################################
+### Compute contributions by simulations (brute-force approach)
+# compute and sort the objective
+Psi = M %*% a;
+Q_sim = quantile( Psi, (1 - c) );
+
+e = mean( abs( a )) / 100; # perturbation
+DQ_simul = matrix( NaN, 1, N) ;
+for( n in 1 : N )
+{
+ # compute gradient
+ a_e = a;
+ a_e[ n ] = a[ n ] + e;
+
+ Psi_e = M %*% a_e;
+ Q_sim_e = quantile(Psi_e, (1 - c) );
+ DQ_simul[ n ] = ( Q_sim_e - Q_sim )/e;
+}
+# compute contributions
+ContrQ_simul = a * t( DQ_simul );
+
+###################################################################################################################
+### Compute contributions analytically
+# compute quantile of standardized marginal (1-dim generator) in simulations this does NOT depend on the allocation...
+gc = quantile(X[ ,1 ], (1 - c));
+
+# ...the dependence on the allocation is analytical
+Q_an = t(Mu) %*% a + gc * sqrt( t(a) %*% Sigma %*% a );
+DQ_an = Mu + gc * Sigma %*% a / sqrt( t(a) %*% Sigma %*% a )[1];
+ContrQ_an = a * DQ_an;
+
+###################################################################################################################
+# plots
+dev.new();
+par( mfrow = c(2,1) );
+bar = barplot( t(ContrQ_an), xlim = c(0, N+1), main = "contributions to VaR, analytical" );
+axis( 1 , at = bar, labels= 1:(N+1));
+
+bar = barplot( t(ContrQ_simul), xlim = c(0, N+1), main = "contributions to VaR, simulations" );
+axis( 1 , at = bar, labels= 1:(N+1));
Added: pkg/Meucci/man/GenerateUniformDrawsOnUnitSphere.Rd
===================================================================
--- pkg/Meucci/man/GenerateUniformDrawsOnUnitSphere.Rd (rev 0)
+++ pkg/Meucci/man/GenerateUniformDrawsOnUnitSphere.Rd 2013-07-16 15:31:33 UTC (rev 2582)
@@ -0,0 +1,33 @@
+\name{GenerateUniformDrawsOnUnitSphere}
+\alias{GenerateUniformDrawsOnUnitSphere}
+\title{Generate a uniform sample on the unit hypersphere, as described in A. Meucci,
+ "Risk and Asset Allocation", Springer, 2005.}
+\usage{
+ GenerateUniformDrawsOnUnitSphere(J, N)
+}
+\arguments{
+ \item{J}{: [scalar] number of draws}
+
+ \item{N}{: [scalar] dimension}
+}
+\value{
+ X : [matrix] (T x N) of draws
+}
+\description{
+ Generate a uniform sample on the unit hypersphere, as
+ described in A. Meucci, "Risk and Asset Allocation",
+ Springer, 2005.
+}
+\note{
+ \item{ Initial script by Xiaoyu Wang - Dec 2006} \item{
+ We decompose X=U*R, where U is a uniform distribution on
+ unit sphere and
+}
+\author{
+ Xavier Valls \email{flamejat at gmail.com}
+}
+\references{
+ \url{http://symmys.com/node/170} See Meucci's script for
+ "GenerateUniformDrawsOnUnitSphere.m"
+}
+
More information about the Returnanalytics-commits
mailing list