[Returnanalytics-commits] r2611 - in pkg/Meucci: . R demo man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sun Jul 21 18:34:54 CEST 2013
Author: xavierv
Date: 2013-07-21 18:34:53 +0200 (Sun, 21 Jul 2013)
New Revision: 2611
Added:
pkg/Meucci/R/PlotMarginalsNormalInverseWishart.R
pkg/Meucci/R/RandNormalInverseWishart.R
pkg/Meucci/demo/S_AnalyzeNormalInverseWishart.R
pkg/Meucci/demo/S_CorrelationPriorUniform.R
pkg/Meucci/demo/S_EvaluationGeneric.R
pkg/Meucci/demo/S_MarkovChainMonteCarlo.R
pkg/Meucci/man/PlotMarginalsNormalInverseWishart.Rd
pkg/Meucci/man/RandNormalInverseWishart.Rd
Modified:
pkg/Meucci/DESCRIPTION
pkg/Meucci/NAMESPACE
Log:
- added all the scripts from chapters 7 and 8 and its associated functions
Modified: pkg/Meucci/DESCRIPTION
===================================================================
--- pkg/Meucci/DESCRIPTION 2013-07-21 10:18:48 UTC (rev 2610)
+++ pkg/Meucci/DESCRIPTION 2013-07-21 16:34:53 UTC (rev 2611)
@@ -30,10 +30,11 @@
xts (>= 0.8),
matlab,
pracma,
- R.utils
+ R.utils,
+ mvtnorm,
+ dlm
Suggests:
quadprog,
- mvtnorm,
limSolve,
Matrix,
MASS,
@@ -83,3 +84,5 @@
'FitExpectationMaximization.R'
'QuantileMixture.R'
'GenerateUniformDrawsOnUnitSphere.R'
+ 'PlotMarginalsNormalInverseWishart.R'
+ 'RandNormalInverseWishart.R'
Modified: pkg/Meucci/NAMESPACE
===================================================================
--- pkg/Meucci/NAMESPACE 2013-07-21 10:18:48 UTC (rev 2610)
+++ pkg/Meucci/NAMESPACE 2013-07-21 16:34:53 UTC (rev 2611)
@@ -28,8 +28,10 @@
export(PartialConfidencePosterior)
export(PerformIidAnalysis)
export(PlotDistributions)
+export(PlotMarginalsNormalInverseWishart)
export(ProjectionStudentT)
export(QuantileMixture)
+export(RandNormalInverseWishart)
export(Raw2Central)
export(Raw2Cumul)
export(RejectOutlier)
Added: pkg/Meucci/R/PlotMarginalsNormalInverseWishart.R
===================================================================
--- pkg/Meucci/R/PlotMarginalsNormalInverseWishart.R (rev 0)
+++ pkg/Meucci/R/PlotMarginalsNormalInverseWishart.R 2013-07-21 16:34:53 UTC (rev 2611)
@@ -0,0 +1,77 @@
+#' Plot the marginals of the normal-inverse-Whishart model.
+#' Described in A. Meucci "Risk and Asset Allocation", Springer, 2005
+#'
+#' @param Mu_Simul : []
+#' @param InvSigma_Simul : []
+#' @param Mu_0 : []
+#' @param T_0 : []
+#' @param Sigma_0 : []
+#' @param Nu_0 : []
+#' @param Legend : []
+#'
+#' @note Numerically and analytically the marginal pdf of
+#' - the first entry of the random vector Mu
+#' - the (1,1)-entry of the random matrix inv(Sigma)
+#' when Mu and Sigma are jointly normal-inverse-Wishart: Mu ~ St(Mu_0,Sigma/T_0)
+#' inv(Sigma) ~ W(Nu_0,inv(Sigma_0)/Nu_0)
+#'
+#' @references
+#' \url{http://symmys.com/node/170}
+#' See Meucci's script for "QuantileMixture.m"
+#'
+#' @author Xavier Valls \email{flamejat@@gmail.com}
+#' @export
+
+PlotMarginalsNormalInverseWishart = function(Mu_Simul, InvSigma_Simul, Mu_0, T_0, Sigma_0, Nu_0, Legend)
+{
+ NumSimulations = nrow( Mu_Simul );
+ NumBins = round( 10 * log( NumSimulations ));
+
+ dev.new();
+
+ #################################################################################################################
+ ### Mu
+ # plot empirical pdf (histogram)
+ par( mfrow = c(2,1) );
+ h = hist(Mu_Simul[ , 1 ], NumBins, plot= FALSE);
+ D = h$mids[ 2 ] - h$mids[ 1 ];
+ n= h$counts /( D * NumSimulations );
+ plot( h$mids, n, type = "h", main = bquote(paste( .(Legend)," ", mu)) );
+ #barplot( n );
+
+ # superimpose analytical expectation
+ points( Mu_0[ 1 ], 0, pch = 21, bg = "red" );
+
+
+ # superimpose analytical pdf
+ x_lo = min(Mu_Simul[ ,1 ]);
+ x_hi = max(Mu_Simul[ ,1 ]);
+ x_grid = seq( x_lo, x_hi, (x_hi-x_lo)/100 );
+ m = Mu_0[ 1 ];
+ s = sqrt(Sigma_0[ 1, 1] / T_0 );
+ f = 1 / s * dt( (x_grid - m )/s, Nu_0 );
+ lines(x_grid, f ,col = "red" );
+
+ #################################################################################################################
+ ### Sigma
+ # plot empirical pdf (histogram)
+ h = hist(InvSigma_Simul[ ,1 ], NumBins, plot= FALSE );
+ D = h$mids[ 2 ] - h$mids[ 1 ];
+ n= h$counts /( D * NumSimulations );
+ plot( h$mids, n, type = "h", main = bquote(paste( .(Legend), " inv(Sigma)")) );
+
+ # superimpose analytical expectation
+ InvSigma_0=solve(Sigma_0);
+ points(InvSigma_0[ 1, 1 ],0, pch = 21, bg = "red" );
+
+
+ # superimpose analytical pdf
+ x_lo = min(InvSigma_Simul[ ,1 ]);
+ x_hi = max(InvSigma_Simul[ ,1 ]);
+ x_grid = seq( x_lo, x_hi, (x_hi-x_lo)/100 );
+ sigma_square = InvSigma_0[ 1, 1] / Nu_0;
+ A = Nu_0 / 2;
+ B = 2 * sigma_square;
+ f = dgamma(x_grid, shape = A, scale = B);
+ lines(x_grid, f, col = "red" );
+}
\ No newline at end of file
Added: pkg/Meucci/R/RandNormalInverseWishart.R
===================================================================
--- pkg/Meucci/R/RandNormalInverseWishart.R (rev 0)
+++ pkg/Meucci/R/RandNormalInverseWishart.R 2013-07-21 16:34:53 UTC (rev 2611)
@@ -0,0 +1,56 @@
+
+#' Generates a multivariate i.i.d. sample of lenght J from the normal-inverse-Wishart distribution, as described in
+#' A. Meucci "Risk and Asset Allocation", Springer, 2005.
+#'
+#' @param Mu_0 : [vector]
+#' @param T_0 : [scalar]
+#' @param Sigma_0 : [matrix]
+#' @param nu_0 : [scalar]
+#' @param J : [scalar]
+#'
+#' @return Mu : [vector]
+#' @return Sigma : [matrix]
+#' @return InvSigma : [matrix]
+#'
+#' @note
+#' Mu|Sigma ~ N(Mu_0,Sigma/T_0)
+#' inv(Sigma) ~ W(Nu_0,inv(Sigma_0)/Nu_0)
+#'
+#' @references
+#' \url{http://symmys.com/node/170}
+#' See Meucci's script for "RandNormalInverseWishart.m"
+#'
+#' @author Xavier Valls \email{flamejat@@gmail.com}
+#' @export
+
+RandNormalInverseWishart = function(Mu_0, T_0, Sigma_0, nu_0, J)
+{
+ N = length( Mu_0 );
+ VecIndex = NULL;
+ for( n in 1 : N )
+ {
+ VecIndex[ n ] = cbind( VecIndex, ( n-1 ) * N +( n:N ) ); ##ok<AGROW>
+ }
+
+ invSigma_0 = solve(Sigma_0) %*% diag( 1, dim( Sigma_0 ));
+ Phi = invSigma_0 / nu_0;
+
+ Mu = NULL;
+ Sigma = NULL;
+ InvSigma = NULL;
+
+ for( j in 1 : J )
+ {
+ Inv_Sigma = rwishart( df = nu_0, Sigma = Phi );
+ InvSigma = rbind( InvSigma, Inv_Sigma[ VecIndex ] );
+
+ S = solve(Inv_Sigma) %*% diag( 1, dim( Inv_Sigma ) );
+ Sigma = rbind( Sigma, S[VecIndex] );
+
+ M = rmvnorm( nrow(Mu_0) * nrow(S/T_0), Mu_0, S/T_0);
+ Mu = rbind( Mu, M );
+ }
+
+ return( list( Mu = Mu, Sigma = Sigma , InvSigma = InvSigma ) )
+}
+
Added: pkg/Meucci/demo/S_AnalyzeNormalInverseWishart.R
===================================================================
--- pkg/Meucci/demo/S_AnalyzeNormalInverseWishart.R (rev 0)
+++ pkg/Meucci/demo/S_AnalyzeNormalInverseWishart.R 2013-07-21 16:34:53 UTC (rev 2611)
@@ -0,0 +1,86 @@
+#' This script familiarizes the users with multivariate Bayesian estimation.
+#' A normal time series is generated a normal-inverse-Wishart prior is set.
+#' The ensuing normal-inverse-Wishart posterior is computed and analyzed numerically and analytically.
+#' Described in A. Meucci,"Risk and Asset Allocation",Springer, 2005, Chapter 7.
+#'
+#' @references
+#' \url{http://symmys.com/node/170}
+#' See Meucci's script for "S_AnalyzeNormalInverseWishart.m"
+#
+#' @author Xavier Valls \email{flamejat@@gmail.com}
+
+###################################################################################################################
+### Inputs
+N = 1; # market dimension
+nSim = 2000;
+
+###################################################################################################################
+### History : X_t ~ N(M,S), t=1,...,T
+# set parameters
+M = 1 * array( 1, N );
+s = 1 * array( 1, N );
+r = 0.7;
+C = (1 - r) * diag( 1, N ) + r * matrix( 1, N, N);
+S = diag( s, N ) %*% C %*% diag( s, N );
+T = 520;
+
+# generate time series
+X = rmvnorm( T, M, S);
+Mu_Hat = apply(X, 2, mean);
+Sigma_Hat = cov(X);
+
+###################################################################################################################
+### Prior : Mu|Sigma ~ N(Mu_0,Sigma/T_0)
+### Omega == inv(Sigma) ~ W(Nu_0,inv(Sigma_0)/Nu_0)
+# set parameters
+Mu_0 = 2 * array( 1, N );
+T_0 = 520;
+s_0 = 2 * array( 1, N );
+r = 0.3;
+C = ( 1 - r ) * diag( 1, N ) + r * matrix( 1, N, N );
+Sigma_0 = diag( s_0, N ) %*% C %*% diag( s_0, N );
+Nu_0 = 520;
+
+# generate simulations
+RNIWPrior= RandNormalInverseWishart(Mu_0, T_0, Sigma_0, Nu_0, nSim);
+
+# plot results
+PlotMarginalsNormalInverseWishart( RNIWPrior$Mu , RNIWPrior$InvSigma, Mu_0, T_0, Sigma_0, Nu_0, "prior" );
+
+###################################################################################################################
+### Posterior : Mu|Sigma ~ N(Mu_1,Sigma/T_1)
+### Omega == inv(Sigma) ~ W(Nu_1,inv(Sigma_1)/Nu_1)
+
+# set parameters
+T_1 = T_0 + T;
+Mu_1 = ( T_0 %*% Mu_0 + T %*% Mu_Hat) / T_1;
+Nu_1 = Nu_0 + T;
+Sigma_1 = ( Nu_0 %*% Sigma_0 + T %*% Sigma_Hat + ( T %*% T_0 / (T + T_0)) %*% (Mu_0 - Mu_Hat) %*% t(Mu_0 - Mu_Hat) ) / Nu_1;
+
+# generate simulations
+RNIWPost = RandNormalInverseWishart(Mu_1, T_1, Sigma_1, Nu_1, nSim);
+
+# plot results
+PlotMarginalsNormalInverseWishart( RNIWPost$Mu, RNIWPost$InvSigma, Mu_1, T_1, Sigma_1, Nu_1, "posterior" );
+
+###################################################################################################################
+### Compute statistics
+Mu_CE_Num = apply( RNIWPost$Mu, 2, mean);
+Mu_CE_Anal = t( Mu_1 );
+Mu_Hat = t( Mu_Hat );
+Mu_0 = t( Mu_0 );
+
+Mu_Scatter_Num = cov( RNIWPost$Mu );
+Mu_Scatter_Anal = Nu_1 / ( Nu_1 - 2 ) * Sigma_1 / T_1;
+
+Sigma_CE_Num = apply(RNIWPost$Sigma,2, mean);
+Sigma_CE_Anal = Sigma_1;
+print(Sigma_Hat);
+print(Sigma_0);
+
+Sigma_Scatter_Num = cov(RNIWPost$Sigma);
+
+InvSigma_CE_Num = apply(RNIWPost$InvSigma, 2, mean);
+S = solve( Sigma_1 );
+InvSigma_CE_Anal = S;
+
Added: pkg/Meucci/demo/S_CorrelationPriorUniform.R
===================================================================
--- pkg/Meucci/demo/S_CorrelationPriorUniform.R (rev 0)
+++ pkg/Meucci/demo/S_CorrelationPriorUniform.R 2013-07-21 16:34:53 UTC (rev 2611)
@@ -0,0 +1,78 @@
+#' This script shows how a jointly uniform prior on the correlations implies that the marginal distribution of
+#' each correlation is peaked around zero , as described in A. Meucci,"Risk and Asset Allocation",Springer, 2005,
+#' Chapter 7.
+#'
+#' @references
+#' \url{http://symmys.com/node/170}
+#' See Meucci's script for "S_CorrelationPriorUniform.m"
+#
+#' @author Xavier Valls \email{flamejat@@gmail.com}
+
+##################################################################################################################
+### Inputs
+N = 3; # dimensionality of the problem
+K = N * (N - 1) / 2;
+J = 10000; # number of simulations
+
+##################################################################################################################
+### Compute correlations in all scenarios
+CorrsAsTensor = array(0, dim = c(J,N,N) );
+Eigs = NULL;
+j = 1;
+
+while( j < J )
+{
+ C = 2 * matrix( runif(K), 1, K ) - 1;
+ Corr = diag( 1, N );
+ k = 0;
+ for( n in 1 : ( N - 1 ) )
+ {
+ for( m in ( n + 1 ) : N )
+ {
+ k = k + 1;
+ Corr[ n, m ] = C[ k ];
+ Corr[ m, n ] = Corr[ n, m ];
+ }
+ }
+ E = eigen(Corr)$values;
+ if( min(E) > 0 )
+ {
+ CorrsAsTensor[ j, , ] = Corr;
+ j = j + 1;
+ }
+ Eigs = rbind(Eigs, t(E) );
+}
+
+#####################################################################################################################
+### Reassemble results in an entry-wise structure that runs on the upper triangular portion of the correlation
+CorrsAsEntries = NULL;
+for( n in 1 : (N-1) )
+{
+ for( m in (n + 1) : N )
+ {
+ el = list( Values = CorrsAsTensor[ , n, m ], Names = paste("{", n, m,"}"));
+ # CorrsAsEntries[ k ]$Names = ["\theta_{") num2str(n) num2str(m) "}"];
+ CorrsAsEntries = rbind( CorrsAsEntries, el )
+ }
+}
+#####################################################################################################################
+### Plots
+# univariate marginals
+K = nrow( CorrsAsEntries );
+Nbins = round( 5 * log( J ) );
+for( k in 1 : K )
+{
+ dev.new();
+ hist(CorrsAsEntries[ k, ]$Values, Nbins, xlab = "", ylab = "", main = bquote( paste("histogram of ", theta, .(CorrsAsEntries[k,]$Names))));
+}
+
+# bivariate marginals
+for( k in 1 : (K-1) )
+{
+ for( j in (k + 1) : K )
+ {
+ dev.new();
+ plot(CorrsAsEntries[k, ]$Values,CorrsAsEntries[j, ]$Values, xlab = "", ylab = "",
+ main = paste( CorrsAsEntries[ k ]$Names, ' - ', CorrsAsEntries[ j ]$Names ));
+ }
+}
Added: pkg/Meucci/demo/S_EvaluationGeneric.R
===================================================================
--- pkg/Meucci/demo/S_EvaluationGeneric.R (rev 0)
+++ pkg/Meucci/demo/S_EvaluationGeneric.R 2013-07-21 16:34:53 UTC (rev 2611)
@@ -0,0 +1,217 @@
+#' Determine the optimal allocation, as described in A. Meucci "Risk and Asset Allocation", Springer, 2005
+#'
+#' @param Market : [struct] market parameters
+#' @param InvestorProfile : [struct] investor's parameters
+#'
+#' @return Allocation : [vector] (N x 1)
+#'
+#' @note
+#' compute optimal allocation, only possible if hidden parameters were known: thus it is not a "decision", we call it a "choice"
+#'
+#' @references
+#' \url{http://symmys.com/node/170}
+#' See Meucci's script for " EvaluationChoiceOptimal.m"
+#'
+#' @author Xavier Valls \email{flamejat@@gmail.com}
+
+
+EvaluationChoiceOptimal = function( Market, InvestorProfile )
+{
+ Exp_Prices = diag( Market$CurrentPrices, length(Market$CurrentPrices) ) %*% ( 1 + Market$LinRets_EV );
+ Cov_Prices = diag( Market$CurrentPrices, length(Market$CurrentPrices) ) %*% Market$LinRets_Cov %*% diag( Market$CurrentPrices, length(Market$CurrentPrices) );
+
+ S = solve( Cov_Prices ) %*% diag( 1, dim(Cov_Prices) );
+ A = (t( Market$CurrentPrices ) %*% S %*% Market$CurrentPrices)[ 1 ];
+ B = (t( Market$CurrentPrices ) %*% S %*% Exp_Prices)[1];
+
+ Gamma = (( InvestorProfile$Budget - InvestorProfile$RiskPropensity * B) / A )[1];
+ Allocation = InvestorProfile$RiskPropensity * S %*% Exp_Prices + Gamma[ 1 ] * S %*% Market$CurrentPrices;
+
+ return( Allocation );
+}
+
+#' Compute the certainty-equivalent statisfaction index , as described in A. Meucci "Risk and Asset Allocation",
+#' Springer, 2005.
+#'
+#' @param Allocation : [vector] (N x 1)
+#' @param Market : [struct] market parameters
+#' @param InvestorProfile : [struct] investor s parameters
+#'
+#' @return CertaintyEquivalent : [scalar]
+#'
+#' @references
+#' \url{http://symmys.com/node/170}
+#' See Meucci's script for " EvaluationSatisfaction.m"
+#'
+#' @author Xavier Valls \email{flamejat@@gmail.com}
+
+
+EvaluationSatisfaction = function( Allocation, Market, InvestorProfile )
+{
+ CertaintyEquivalent = t(Allocation) %*% diag( Market$CurrentPrices, length( Market$CurrentPrices ) ) %*% (1 + Market$LinRets_EV) - 1 / (2 * InvestorProfile$RiskPropensity) * t( Allocation ) %*% diag( Market$CurrentPrices, length( Market$CurrentPrices )) %*% Market$LinRets_Cov %*% diag( Market$CurrentPrices, length( Market$CurrentPrices )) %*% Allocation ;
+
+ return( CertaintyEquivalent[1] )
+}
+
+
+#' Determine the allocation of the best performer, as described in A. Meucci "Risk and Asset Allocation",
+#' Springer, 2005.
+#'
+#' @param Market : [struct] market parameters
+#' @param InvestorProfile : [struct] investors parameters
+#'
+#' @return Allocation : [vector] (N x 1)
+#'
+#' @note
+#' scenario-dependent decision that tries to pick the optimal allocation
+#'
+#' @references
+#' \url{http://symmys.com/node/170}
+#' See Meucci's script for "EvaluationDecisionBestPerformer.m"
+#'
+#' @author Xavier Valls \email{flamejat@@gmail.com}
+
+EvaluationDecisionBestPerformer = function( Market, InvestorProfile )
+{
+ # find index of best performer
+ B = which.max( Market$LinRetsSeries[ nrow(Market$LinRetsSeries) , ] ); ##ok<ASGLU>
+
+ # invest in that asset
+ I = diag( 1, length(Market$CurrentPrices) );
+ Allocation = InvestorProfile$Budget * I[ , B ] / Market$CurrentPrices[ B ];
+
+ return( Allocation );
+}
+
+
+#' Determine the cost of allocation, as described in A. Meucci "Risk and Asset Allocation", Springer, 2005.
+#'
+#' @param Allocation : [vector] (N x 1)
+#' @param Market : [struct] market parameters
+#' @param InvestorProfile : [struct] investor's parameters
+#'
+#' @return C_Plus : [scalar] cost
+#'
+#' @references
+#' \url{http://symmys.com/node/170}
+#' See Meucci's script for "EvaluationDecisionBestPerformer.m"
+#'
+#' @author Xavier Valls \email{flamejat@@gmail.com}
+
+EvaluationCost = function( Allocation, Market, InvestorProfile )
+{
+ aXi = t(Allocation) %*% diag( Market$CurrentPrices, length( Market$CurrentPrices ) ) %*% (1 + Market$LinRets_EV);
+ aPhia = t(Allocation) %*% diag( Market$CurrentPrices, length( Market$CurrentPrices ) ) %*% Market$LinRets_Cov %*% diag( Market$CurrentPrices, length( Market$CurrentPrices ) ) %*% Allocation;
+
+ C = ( 1 - InvestorProfile$BaR ) * InvestorProfile$Budget - aXi + sqrt(2 %*% aPhia) * erfinv( 2 * InvestorProfile$Confidence - 1);
+ C_Plus = max(C, 0);
+ return( C_Plus );
+
+}
+
+
+#' This script evaluates a generic allocation decision (in this case the "best performer" strategy, which fully
+#' invest the budget in the last period's best performer).
+#' It displays the distribution of satisfaction, cost of constraint violation and opportunity cost for each value
+#' of the market stress-test parameters (in this case the correlation).
+#' Described in A. Meucci "Risk and Asset Allocation", Springer, 2005, Chapter 8.
+#'
+#' @references
+#' \url{http://symmys.com/node/170}
+#' See Meucci's script for "S_EvaluationGeneric.m"
+#
+#' @author Xavier Valls \email{flamejat@@gmail.com}
+
+##################################################################################################################
+### Inputs
+NumScenarios = 1000;
+NumCorrelations = 5;
+Bottom_Correlation = 0;
+Top_Correlation = 0.99;
+
+##################################################################################################################
+### Input investor's parameters
+InvestorProfile = NULL;
+InvestorProfile$Budget = 10000;
+InvestorProfile$RiskPropensity = 30;
+InvestorProfile$Confidence = 0.9;
+InvestorProfile$BaR = 0.2;
+
+##################################################################################################################
+### Input market parameters
+NumAssets = 10;
+a = 0.5; # effect of correlation on expected values and volatility (hidden)
+Bottom = 0.06;
+Top = 0.36;
+Step = (Top - Bottom) / (NumAssets - 1);
+v = seq( Bottom, Top, Step ) ; # volatility vector
+Market = NULL;
+Market$T = 20; # not hidden
+Market$CurrentPrices = 10 * array( 1, NumAssets); # not hidden
+
+##################################################################################################################
+Step = (Top_Correlation - Bottom_Correlation) / (NumCorrelations - 1);
+Overall_Correlations = seq( Bottom_Correlation, Top_Correlation, Step );
+
+Suboptimal = NULL;
+Suboptimal$StrsTst_Satisfaction = NULL;
+Suboptimal$StrsTst_CostConstraints = NULL;
+Suboptimal$StrsTst_OppCost = NULL;
+Optimal = NULL;
+Optimal$StrsTst_Satisfaction = NULL;
+
+for( t in 1 : length(Overall_Correlations) )
+{
+ # input the (hidden) market parameters (only correlations, we assume standard deviations and expected values fixed and known)
+ Market$St_Devations = ( 1 + a * Overall_Correlations[ t ]) * v; # hidden
+ Market$LinRets_EV = 0.5 * Market$St_Devations; # hidden
+
+ Correlation = ( 1 - Overall_Correlations[ t ] ) * diag( 1, NumAssets) + Overall_Correlations[ t ] * matrix( 1, NumAssets, NumAssets);
+ Market$LinRets_Cov = diag( Market$St_Devations, length(Market$St_Devations) ) %*% Correlation %*% diag( Market$St_Devations, length(Market$St_Devations) )
+
+
+ ##################################################################################################################
+ # compute optimal allocation, only possible if hidden parameters were known: thus it is not a "decision", we call it a "choice"
+ Allocation = EvaluationChoiceOptimal( Market, InvestorProfile );
+ Satisfaction_Optimal = EvaluationSatisfaction( Allocation, Market, InvestorProfile );
+
+ ##################################################################################################################
+ # choose allocation based on available information
+ StrsTst_TrueSatisfaction = NULL;
+ StrsTst_CostConstraints = NULL;
+
+ for( s in 1 : NumScenarios )
+ {
+ # generate scenarios i_T of information I_T
+ Market$LinRetsSeries = rmvnorm( Market$T, Market$LinRets_EV, Market$LinRets_Cov );
+
+ # scenario-dependent decision that tries to pick the optimal allocation
+ Allocation = EvaluationDecisionBestPerformer( Market, InvestorProfile );
+ TrueSatisfaction = EvaluationSatisfaction( Allocation, Market, InvestorProfile );
+ CostConstraints = EvaluationCost( Allocation, Market, InvestorProfile );
+
+ StrsTst_TrueSatisfaction = cbind( StrsTst_TrueSatisfaction, TrueSatisfaction ); ##ok<*AGROW>
+ StrsTst_CostConstraints = cbind( StrsTst_CostConstraints, CostConstraints );
+ }
+
+ Suboptimal$StrsTst_CostConstraints = rbind( Suboptimal$StrsTst_CostConstraints, StrsTst_CostConstraints );
+ Suboptimal$StrsTst_Satisfaction = rbind( Suboptimal$StrsTst_Satisfaction, StrsTst_TrueSatisfaction );
+ Suboptimal$StrsTst_OppCost = rbind( Suboptimal$StrsTst_OppCost, Satisfaction_Optimal - StrsTst_TrueSatisfaction + StrsTst_CostConstraints );
+ Optimal$StrsTst_Satisfaction = rbind( Optimal$StrsTst_Satisfaction, Satisfaction_Optimal );
+}
+
+##################################################################################################################
+### Display
+NumVBins = round(10 * log(NumScenarios));
+
+# optimal allocation vs. allocation decision
+for( t in 1 : length(Overall_Correlations) )
+{
+ dev.new();
+ par( mfrow = c( 3, 1) )
+ hist(Suboptimal$StrsTst_Satisfaction[ t, ], NumVBins, main = "satisfaction", xlab ="", ylab = "" );
+
+ hist(Suboptimal$StrsTst_CostConstraints[ t, ], NumVBins, main = "constraint violation cost", xlab ="", ylab = "");
+
+ hist(Suboptimal$StrsTst_OppCost[ t, ], NumVBins, main = "opportunity cost", xlab ="", ylab = "");
+}
\ No newline at end of file
Added: pkg/Meucci/demo/S_MarkovChainMonteCarlo.R
===================================================================
--- pkg/Meucci/demo/S_MarkovChainMonteCarlo.R (rev 0)
+++ pkg/Meucci/demo/S_MarkovChainMonteCarlo.R 2013-07-21 16:34:53 UTC (rev 2611)
@@ -0,0 +1,52 @@
+#' This script illustrates the Metropolis-Hastings algorithm, as described in A. Meucci,"Risk and Asset Allocation",
+#' Springer, 2005, Chapter 7.
+#'
+#' @references
+#' \url{http://symmys.com/node/170}
+#' See Meucci's script for "S_MarkovChainMonteCarlo.m"
+#
+#' @author Xavier Valls \email{flamejat@@gmail.com}
+
+##################################################################################################################
+### Set-up target and candidate
+# kernel of the target distribution
+kernel = function(x) dnorm( x, 1, 0.5 );
+
+# parameters of the normal candidate distribution
+mu = 0;
+sig = 5;
+
+##################################################################################################################
+### Set up MH algorithm
+nSim = 10000;
+xt = matrix( NaN, nSim, 1);
+xt[ 1 ] = 0;
+nacc = 0;
+for( i in 2 : nSim )
+{
+ # normal candidate
+ r = mu + sig * rnorm(1);
+ # kernel at candidate
+ f1 = kernel( r );
+ # kernel at past
+ f2 = kernel( xt[ i-1 ] );
+ prob = f1 / f2;
+ xt[ i ] = xt[ i-1 ];
+ if( prob > 1 || runif(1) > (1 - prob) )
+ {
+ xt[ i ] = r;
+ nacc = nacc + 1;
+ }
+}
+##################################################################################################################
+### Post-process output
+# acceptance rate
+print( nacc / nSim );
+
+# display MCMC over time
+dev.new();
+plot( xt, type = "l" );
+
+# distribution
+dev.new();
+hist( xt, log(10*nSim) );
Added: pkg/Meucci/man/PlotMarginalsNormalInverseWishart.Rd
===================================================================
--- pkg/Meucci/man/PlotMarginalsNormalInverseWishart.Rd (rev 0)
+++ pkg/Meucci/man/PlotMarginalsNormalInverseWishart.Rd 2013-07-21 16:34:53 UTC (rev 2611)
@@ -0,0 +1,43 @@
+\name{PlotMarginalsNormalInverseWishart}
+\alias{PlotMarginalsNormalInverseWishart}
+\title{Plot the marginals of the normal-inverse-Whishart model.
+Described in A. Meucci "Risk and Asset Allocation", Springer, 2005}
+\usage{
+ PlotMarginalsNormalInverseWishart(Mu_Simul,
+ InvSigma_Simul, Mu_0, T_0, Sigma_0, Nu_0, Legend)
+}
+\arguments{
+ \item{Mu_Simul}{: []}
+
+ \item{InvSigma_Simul}{: []}
+
+ \item{Mu_0}{: []}
+
+ \item{T_0}{: []}
+
+ \item{Sigma_0}{: []}
+
+ \item{Nu_0}{: []}
+
+ \item{Legend}{: []}
+}
+\description{
+ Plot the marginals of the normal-inverse-Whishart model.
+ Described in A. Meucci "Risk and Asset Allocation",
+ Springer, 2005
+}
+\note{
+ Numerically and analytically the marginal pdf of - the
+ first entry of the random vector Mu - the (1,1)-entry of
+ the random matrix inv(Sigma) when Mu and Sigma are
+ jointly normal-inverse-Wishart: Mu ~ St(Mu_0,Sigma/T_0)
+ inv(Sigma) ~ W(Nu_0,inv(Sigma_0)/Nu_0)
+}
+\author{
+ Xavier Valls \email{flamejat at gmail.com}
+}
+\references{
+ \url{http://symmys.com/node/170} See Meucci's script for
+ "QuantileMixture.m"
+}
+
Added: pkg/Meucci/man/RandNormalInverseWishart.Rd
===================================================================
--- pkg/Meucci/man/RandNormalInverseWishart.Rd (rev 0)
+++ pkg/Meucci/man/RandNormalInverseWishart.Rd 2013-07-21 16:34:53 UTC (rev 2611)
@@ -0,0 +1,42 @@
+\name{RandNormalInverseWishart}
+\alias{RandNormalInverseWishart}
+\title{Generates a multivariate i.i.d. sample of lenght J from the normal-inverse-Wishart distribution, as described in
+A. Meucci "Risk and Asset Allocation", Springer, 2005.}
+\usage{
+ RandNormalInverseWishart(Mu_0, T_0, Sigma_0, nu_0, J)
+}
+\arguments{
+ \item{Mu_0}{: [vector]}
+
+ \item{T_0}{: [scalar]}
+
+ \item{Sigma_0}{: [matrix]}
+
+ \item{nu_0}{: [scalar]}
+
+ \item{J}{: [scalar]}
+}
+\value{
+ Mu : [vector]
+
+ Sigma : [matrix]
+
+ InvSigma : [matrix]
+}
+\description{
+ Generates a multivariate i.i.d. sample of lenght J from
+ the normal-inverse-Wishart distribution, as described in
+ A. Meucci "Risk and Asset Allocation", Springer, 2005.
+}
+\note{
+ Mu|Sigma ~ N(Mu_0,Sigma/T_0) inv(Sigma) ~
+ W(Nu_0,inv(Sigma_0)/Nu_0)
+}
+\author{
+ Xavier Valls \email{flamejat at gmail.com}
+}
+\references{
+ \url{http://symmys.com/node/170} See Meucci's script for
+ "RandNormalInverseWishart.m"
+}
+
More information about the Returnanalytics-commits
mailing list