[Returnanalytics-commits] r2571 - pkg/Meucci/demo
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sun Jul 14 22:27:53 CEST 2013
Author: xavierv
Date: 2013-07-14 22:27:52 +0200 (Sun, 14 Jul 2013)
New Revision: 2571
Added:
pkg/Meucci/demo/S_EstimateQuantileEvaluation.R
pkg/Meucci/demo/S_ShrinkageEstimators.R
Log:
-added S_ShrinkageEstimators and S_EstimateQuantileEvaluation demo scripts
Added: pkg/Meucci/demo/S_EstimateQuantileEvaluation.R
===================================================================
--- pkg/Meucci/demo/S_EstimateQuantileEvaluation.R (rev 0)
+++ pkg/Meucci/demo/S_EstimateQuantileEvaluation.R 2013-07-14 20:27:52 UTC (rev 2571)
@@ -0,0 +1,150 @@
+#'This script familiarizes the user with the evaluation of an estimator:replicability, loss, error,
+#'bias and inefficiency as described in A. Meucci,"Risk and Asset Allocation", Springer, 2005, Chapter 4.
+#'
+#' \url{http://symmys.com/node/170}
+#' See Meucci's script for "S_EstimateQuantileEvaluation.m"
+#'
+#' @author Xavier Valls \email{flamejat@@gmail.com}
+
+##################################################################################################################
+### Inputs
+
+T = 52;
+a = 0.5;
+m_Y = 0.1;
+s_Y = 0.2;
+m_Z = 0;
+s_Z = 0.15;
+
+##################################################################################################################
+### plain estimation
+
+# functional of the distribution to be estimated
+G_fX = QuantileMixture( 0.5, a, m_Y, s_Y, m_Z, s_Z);
+print( G_fX );
+
+# series generated by "nature": do not know the distribution
+P = runif( T );
+i_T = QuantileMixture( P, a, m_Y, s_Y, m_Z, s_Z );
+
+G_Hat_e = function(X) apply( X, 1, median );
+G_Hat_b = function(X) apply( X, 1, mean );
+
+Ge = G_Hat_e( t( as.matrix( i_T )) ); # tentative estimator of unknown functional
+Gb = G_Hat_b( t( as.matrix( i_T )) ); # tentative estimator of unknown functional
+print(Ge);
+print(Gb);
+
+##################################################################################################################
+### replicability vs. "luck"
+
+# functional of the distribution to be estimated
+G_fX = QuantileMixture( 0.5, a, m_Y, s_Y, m_Z, s_Z);
+
+# randomize series generated by "nature" to check replicability
+nSim = 10000;
+I_T = c();
+
+for( t in 1:T )
+{
+ P = runif(nSim);
+ Simul = QuantileMixture(P,a,m_Y,s_Y,m_Z,s_Z);
+ I_T = cbind( I_T, Simul );
+}
+
+Ge = G_Hat_e( I_T ); # tentative estimator of unknown functional
+Gb = G_Hat_b( I_T ); # tentative estimator of unknown functional
+
+Loss_Ge = ( Ge - G_fX ) ^ 2;
+Loss_Gb = ( Gb - G_fX ) ^ 2;
+
+Err_Ge = sqrt( mean( Loss_Ge));
+Err_Gb = sqrt( mean( Loss_Gb));
+
+Bias_Ge = abs(mean(Ge)-G_fX);
+Bias_Gb = abs(mean(Gb)-G_fX);
+
+Ineff_Ge = sd(Ge);
+Ineff_Gb = sd(Gb);
+
+###################################################################################################################
+dev.new();
+NumBins = round( 10 * log( nSim ));
+par( mfrow = c( 2, 1) );
+hist( Ge, NumBins, main = "estimator e" );
+points(G_fX, 0, pch = 21, bg = "red" );
+
+hist(Gb, NumBins, main = "estimator b" );
+points(G_fX, 0, pch = 21, bg = "red" );
+
+#loss
+dev.new();
+par( mfrow = c( 2, 1) );
+hist(Loss_Ge, NumBins, main = "loss of estimator e" );
+hist(Loss_Gb, NumBins, main = "loss of estimator b" );
+
+###################################################################################################################
+### stress test replicability
+m_s = seq( 0, 0.2, 0.02 );
+
+Err_Gesq = NULL; Bias_Gesq = NULL; Ineff_Gesq = NULL;
+Err_Gbsq = NULL; Bias_Gbsq = NULL; Ineff_Gbsq = NULL;
+
+for( j in 1 : length( m_s ) )
+{
+ m_Y = m_s[ j ];
+ # functional of the distribution to be estimated
+ G_fX = QuantileMixture( 0.5, a, m_Y, s_Y, m_Z, s_Z );
+
+ # randomize series generated by "nature" to check replicability
+ nSim = 10000;
+ I_T = NULL;
+ for( t in 1 : T )
+ {
+ P = runif(nSim);
+ Simul = QuantileMixture(P, a, m_Y, s_Y, m_Z, s_Z);
+ I_T = cbind( I_T, Simul );
+ }
+
+ Ge = G_Hat_e( I_T ); # tentative estimator of unknown functional
+ Gb = G_Hat_b( I_T ); # tentative estimator of unknown functional
+
+ Loss_Ge = ( Ge - G_fX ) ^ 2;
+ Loss_Gb = ( Gb - G_fX ) ^ 2;
+
+ Err_Ge = sqrt( mean( Loss_Ge ) );
+ Err_Gb = sqrt( mean( Loss_Gb ) );
+
+ Bias_Ge = abs( mean( Ge ) - G_fX );
+ Bias_Gb = abs( mean( Gb ) - G_fX );
+
+ Ineff_Ge = std( Ge );
+ Ineff_Gb = std( Gb );
+
+ #store results
+ Err_Gesq = cbind( Err_Gesq, Err_Ge ^ 2); ##ok<*AGROW>
+ Err_Gbsq = cbind(Err_Gbsq, Err_Gb^2);
+
+ Bias_Gesq = cbind( Bias_Gesq, Bias_Ge^2 );
+ Bias_Gbsq = cbind( Bias_Gbsq, Bias_Gb^2 );
+
+ Ineff_Gesq = cbind( Ineff_Gesq, Ineff_Ge ^ 2 );
+ Ineff_Gbsq = cbind( Ineff_Gbsq, Ineff_Gb ^ 2 );
+}
+
+###################################################################################################################
+### printlay results
+dev.new();
+par( mfrow = c( 2, 1) );
+b = barplot(Bias_Gesq +Ineff_Gesq , col = "red", main = "stress-test of estimator e");
+barplot( Ineff_Gesq, col="blue", add = TRUE);
+lines( b, Err_Gesq);
+legend( "topleft", 1.9, c( "bias²", "ineff²", "error²" ), col = c( "red","blue", "black" ),
+ lty=1, lwd=c(5,5,1),bg = "gray90" );
+
+
+b = barplot(Bias_Gbsq+Ineff_Gbsq, col = "red", main = "stress-test of estimator b");
+barplot( Ineff_Gbsq, col="blue", add = TRUE);
+lines( b, Err_Gbsq);
+legend( "topleft", 1.9, c( "bias²", "ineff²", "error²" ), col = c( "red","blue", "black" ),
+ lty=1, lwd=c(5,5,1),bg = "gray90" );
\ No newline at end of file
Added: pkg/Meucci/demo/S_ShrinkageEstimators.R
===================================================================
--- pkg/Meucci/demo/S_ShrinkageEstimators.R (rev 0)
+++ pkg/Meucci/demo/S_ShrinkageEstimators.R 2013-07-14 20:27:52 UTC (rev 2571)
@@ -0,0 +1,70 @@
+library(mvtnorm);
+
+#' This script computes the multivariate shrinkage estimators of location and scatter under the normal assumption,
+#' as described in A. Meucci, "Risk and Asset Allocation", Springer, 2005, Chapter 4.
+#'
+#' @references
+#' \url{http://symmys.com/node/170}
+#' See Meucci's script for "S_ShrinkageEstimators.m"
+#'
+
+##################################################################################################################
+### Inputs
+N = 5;
+T = 30;
+Mu = runif(N);
+A = matrix(runif( N* N), N, N) - 0.5;
+Sigma = A %*% t(A);
+
+#################################################################################################################
+### Generate normal sample
+X = rmvnorm( T, Mu, Sigma );
+
+# estimate sample parameters
+Mu_hat = matrix( apply( X, 2, mean ));
+Sigma_hat = cov( X ) * (T - 1) / T;
+
+#################################################################################################################
+### Shrinkage of location parameter
+
+# target
+b = matrix( 0, N, 1 );
+
+# compute optimal weight
+Lambda_hat = eigen( Sigma_hat )$values;
+a = 1 / T * ( sum( Lambda_hat ) - 2 * max( Lambda_hat) ) / ( t( Mu_hat-b) %*% ( Mu_hat - b ) );
+
+# restrict to sensible weight
+a = max( 0, min( a, 1) );
+
+# shrink
+Mu_shr = ( 1 - a ) * Mu_hat + a * b;
+
+print(Mu_hat);
+print(Mu_shr);
+
+#################################################################################################################
+### Shrinkage of scatter parameter
+
+# target
+C = mean(Lambda_hat) * diag( 1, N );
+
+# compute optimal weight
+Numerator = 0;
+for( t in 1 : T )
+{
+ Numerator = Numerator + 1 / T * sum(diag( (matrix(X[ t, ])%*%X[ t, ] - Sigma_hat ) ^ 2 ));
+}
+
+Denominator = sum( diag( ( Sigma_hat - C ) ^ 2 ));
+a = 1 / T * Numerator / Denominator;
+
+# restrict to sensible weight
+a = max(0, min(a, 1));
+
+# shrink
+Sigma_shr = ( 1 - a ) * Sigma_hat + a * C;
+
+print(Sigma_hat);
+print(Sigma_shr);
+
More information about the Returnanalytics-commits
mailing list