[Returnanalytics-commits] r2740 - in pkg/Meucci: R demo

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

Author: xavierv
Date: 2013-08-08 12:24:51 +0200 (Thu, 08 Aug 2013)
New Revision: 2740

Added:
   pkg/Meucci/demo/S_MultiVarSqrRootRule.R
Modified:
   pkg/Meucci/R/TwoDimEllipsoid.R
   pkg/Meucci/demo/S_CrossSectionIndustries.R
Log:
- added S_MultiVarSqrRootRule demo script from chapter 3 and fixed function TwoDimEllipsoid

Modified: pkg/Meucci/R/TwoDimEllipsoid.R
===================================================================
--- pkg/Meucci/R/TwoDimEllipsoid.R	2013-08-07 22:47:20 UTC (rev 2739)
+++ pkg/Meucci/R/TwoDimEllipsoid.R	2013-08-08 10:24:51 UTC (rev 2740)
@@ -32,8 +32,8 @@
 	for( i in 1 : NumSteps )
 	{
 	    # normalized variables (parametric representation of the ellipsoid)
-	    y = c( cos( Angle[ i ] ), sin( Angle[ i ] ) );
-	    Centered_Ellipse = c( Centered_Ellipse, Eigen$vectors %*% diag(sqrt(Eigen$values)) %*% y );   ##ok<AGROW>
+	    y = rbind( cos( Angle[ i ] ), sin( Angle[ i ] ) );
+	    Centered_Ellipse = c( Centered_Ellipse, Eigen$vectors %*% diag(sqrt(Eigen$values), length(Eigen$values)) %*% y );   ##ok<AGROW>
 	}
 
 	R = Location %*% array( 1, NumSteps ) + Scale * Centered_Ellipse;

Modified: pkg/Meucci/demo/S_CrossSectionIndustries.R
===================================================================
--- pkg/Meucci/demo/S_CrossSectionIndustries.R	2013-08-07 22:47:20 UTC (rev 2739)
+++ pkg/Meucci/demo/S_CrossSectionIndustries.R	2013-08-08 10:24:51 UTC (rev 2740)
@@ -10,7 +10,7 @@
 ##################################################################################################################
 ### Load data
 # loads weekly stock returns X and indices stock returns F
-load(" ../data/securitiesTS.Rda");
+load("../data/securitiesTS.Rda");
 Data_Securities = securitiesTS$data[ , -1 ]; # 1st column is date
 
 load("../data/securitiesIndustryClassification.Rda");

Added: pkg/Meucci/demo/S_MultiVarSqrRootRule.R
===================================================================
--- pkg/Meucci/demo/S_MultiVarSqrRootRule.R	                        (rev 0)
+++ pkg/Meucci/demo/S_MultiVarSqrRootRule.R	2013-08-08 10:24:51 UTC (rev 2740)
@@ -0,0 +1,69 @@
+#' This script illustrates the multivariate square root rule-of-thumb 
+#' Described in A. Meucci,"Risk and Asset Allocation", Springer, 2005,  Chapter 3.
+#'
+#' @references
+#' \url{http://symmys.com/node/170}
+#' See Meucci's script for "S_MultiVarSqrRootRule.m"
+#
+#' @author Xavier Valls \email{flamejat@@gmail.com}
+
+##################################################################################################################
+### Load data
+load("../data/swaps.Rda");
+
+##################################################################################################################
+### Aggregation steps in days
+Steps = c( 1, 5, 22 ); 
+
+##################################################################################################################
+### Plots
+Agg = list();
+#names(Agg)=c( "M_hat" , "S_hat", "M_norm", "S_norm");
+
+dev.new();
+plot( swaps$X[ , 1 ],swaps$X[ , 2], xlab = swaps$Names[[1]][1], ylab = swaps$Names[[2]][1] );
+T = nrow(swaps$X );
+
+for( s in 1 : length(Steps) )
+{
+
+    
+
+    # compute series at aggregated time steps
+    k = Steps[ s ];
+    AggX = NULL;
+    t = 1;
+    while( ( t + k + 1 ) <= T )
+    {
+        NewTerm = apply( matrix(swaps$X[ t : (t+k-1), ], ,ncol(swaps$X) ),2,sum);
+        AggX = rbind( AggX, NewTerm ); ##ok<AGROW>
+        t = t + k;
+    }
+
+
+    # empirical mean/covariance
+    
+    if(s==1)
+    {
+        M_hat = matrix(apply(AggX, 2, mean));
+        S_hat = cov(AggX);
+        Agg[[s]] = list( M_hat = M_hat, S_hat = S_hat, M_norm = k / Steps[ 1 ] * M_hat,  S_norm = k / Steps[ 1 ] * S_hat );
+    }else
+    {
+
+        Agg[[s]] = list(
+                M_hat = matrix(apply(AggX,2,mean)),
+                S_hat = cov(AggX),
+                # mean/covariance implied by propagation law of risk for invariants
+                M_norm = k / Steps[ 1 ] * Agg[[ 1 ]]$M_hat, 
+                S_norm = k / Steps[ 1 ] * Agg[[ 1 ]]$S_hat
+                );
+    }
+
+    # plots
+    h1 = TwoDimEllipsoid( Agg[[ s ]]$M_norm, Agg[[ s ]]$S_norm, 1, 0, 0 );
+    
+    h2 = TwoDimEllipsoid( Agg[[ s ]]$M_hat, Agg[[ s ]]$S_hat, 1, 0, 0 );
+
+    
+}