[Returnanalytics-commits] r2520 - pkg/Meucci/demo

[noreply at r-forge.r-project.org]

Author: xavierv
Date: 2013-07-08 09:29:30 +0200 (Mon, 08 Jul 2013)
New Revision: 2520

Added:
   pkg/Meucci/demo/S_EigenvalueDispersion.R
Log:
-added S_EigenvalueDispersion demo script

Added: pkg/Meucci/demo/S_EigenvalueDispersion.R
===================================================================
--- pkg/Meucci/demo/S_EigenvalueDispersion.R	                        (rev 0)
+++ pkg/Meucci/demo/S_EigenvalueDispersion.R	2013-07-08 07:29:30 UTC (rev 2520)
@@ -0,0 +1,60 @@
+library(mvtnorm)
+
+#'This script displays the sample eigenvalues dispersion phenomenon, 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_EigenValueDispersion.R"
+#'
+#' @author Xavier Valls \email{flamejat@@gmail.com}
+
+##################################################################################################################
+### Inputs
+
+N = 50;
+SampleLenght = seq( N , 10 * N, N)
+nSim = 50;
+
+##################################################################################################################
+### Generate mesh and surface
+
+Mu = matrix( 0, N, 1 );
+Sigma= diag( 1, N );
+
+# compute true eigenvalues
+Eigen = eigen(Sigma);
+Index = order( -( Eigen$values ));
+EVec  = Eigen$vectors[ , Index ];
+EVal  = diag( Eigen$Values[ Index, Index ]);
+
+# compute eigenvalues of sample estimator
+nSampleLenght = length( SampleLenght );
+Store_EVal_Hat = matrix( NaN, nSampleLenght, N ); # preallocation for speed
+for( i in 1 : nSampleLenght )
+ {  
+    T = SampleLenght[ i ];
+    EVal_Hat = 0;
+    for( n in 1 : nSim )
+    {
+        X = rmvnorm( T, Mu, Sigma );
+        Sigma_Hat = cov( X );
+        L = eigen( Sigma_Hat )$values;
+        Index = order(-(L));
+        L = L[ Index];
+
+        EVal_Hat = EVal_Hat + L;
+    }
+    EVal_Hat = EVal_Hat / nSim;
+
+    Store_EVal_Hat[ i, ] = t(EVal_Hat);
+}
+
+##################################################################################################################
+### Display surface
+dev.new();
+
+persp( SampleLenght/N, 1 :N , Store_EVal_Hat,
+    theta = 7 * 45, phi = 30, expand=0.6, col='lightblue', shade=0.75, ltheta=120, 
+    ticktype='detailed', xlab = "eigenvalue #", ylab = "sample lenght/N");
+