[Returnanalytics-commits] r2749 - in pkg/Meucci: data demo

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

Author: xavierv
Date: 2013-08-09 13:32:24 +0200 (Fri, 09 Aug 2013)
New Revision: 2749

Added:
   pkg/Meucci/data/sectorsTS.Rda
   pkg/Meucci/demo/S_TimeSeriesIndustries.R
   pkg/Meucci/demo/S_TimeSeriesVsCrossSectionIndustries.R
Log:
- added S_TimeseriesIndustries and S_TimeseriesVsCrossSectionIndustries demo scripts from Chapter 3

Added: pkg/Meucci/data/sectorsTS.Rda
===================================================================
(Binary files differ)


Property changes on: pkg/Meucci/data/sectorsTS.Rda
___________________________________________________________________
Added: svn:mime-type
   + application/octet-stream

Added: pkg/Meucci/demo/S_TimeSeriesIndustries.R
===================================================================
--- pkg/Meucci/demo/S_TimeSeriesIndustries.R	                        (rev 0)
+++ pkg/Meucci/demo/S_TimeSeriesIndustries.R	2013-08-09 11:32:24 UTC (rev 2749)
@@ -0,0 +1,76 @@
+#' This script fits a time-series linear factor computing the industry factors loadings, as 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_TimeSeriesIndustries.m"
+#'
+#' @author Xavier Valls \email{flamejat@@gmail.com}
+
+##################################################################################################################
+### Loads weekly stock returns X and indices stock returns F
+load("../data/securitiesTS.Rda");
+Data_Securities = securitiesTS$data[ , -1 ]; # 1st column is date
+
+load("../data/sectorsTS.Rda");
+Data_Sectors = sectorsTS$data[ , -(1:2) ]; #1st column is for date, 2nd column is SPX index
+
+##################################################################################################################
+### Estimation
+# linear returns for stocks
+X = diff( Data_Securities ) / Data_Securities[ -nrow(Data_Securities), ];
+
+# linear return for the factors
+F = diff(Data_Sectors) / Data_Sectors[ -nrow(Data_Sectors), ];
+
+T = dim(X)[1];
+N = dim(X)[2];
+K = dim(F)[ 2 ];
+
+# compute sample estimates
+E_X = matrix( apply(X, 2, mean ) );
+E_F = matrix( apply(F, 2, mean ) );
+
+XF = cbind( X, F )
+SigmaJoint_XF = (dim(XF)[1]-1)/dim(XF)[1] * cov(XF);
+Sigma_X       = SigmaJoint_XF[ 1:N, 1:N ];
+Sigma_XF      = SigmaJoint_XF[ 1:N, (N+1):(N+K) ];
+Sigma_F       = SigmaJoint_XF[ (N+1):(N+K), (N+1):(N+K) ];
+Corr_F  	  = cov2cor(Sigma_F);
+Corr_F  	  = tril(Corr_F, -1);
+
+# compute OLS loadings for the linear return model
+X_ = X - repmat( t(E_X), T, 1 );
+F_ = F - repmat( t(E_F), T, 1 );
+B  = Sigma_XF %*% solve(Sigma_F);
+U  = X_ - F_ %*% t(B);
+
+##################################################################################################################
+### Residual analysis
+
+UF = cbind(U,F);
+SigmaJoint_UF = ( dim( UF )[1]-1 )/dim( UF )[1] * cov( UF );
+CorrJoint_UF  = cov2cor( SigmaJoint_UF );
+
+# correlations of residuals with factors is null
+Corr_UF      = CorrJoint_UF[ 1:N, (N+1):(N+K) ];
+mean_Corr_UF = mean( abs( as.array( Corr_UF )) );
+max_Corr_UF  = max( abs( as.array( Corr_UF )) );
+
+disp( mean_Corr_UF );
+disp( max_Corr_UF );
+
+dev.new();
+hist( Corr_UF, 100);
+
+# correlations between residuals is not null
+Corr_U = tril( CorrJoint_UF[ 1:N, 1:N ], -1 );
+Corr_U = Corr_U[ Corr_U != 0 ];
+mean_Corr_U = mean( abs( Corr_U ) );
+max_Corr_U  = max( abs( Corr_U ) );
+disp( mean_Corr_U );
+disp( max_Corr_U );
+
+dev.new();
+hist(Corr_U, 100);
+

Added: pkg/Meucci/demo/S_TimeSeriesVsCrossSectionIndustries.R
===================================================================
--- pkg/Meucci/demo/S_TimeSeriesVsCrossSectionIndustries.R	                        (rev 0)
+++ pkg/Meucci/demo/S_TimeSeriesVsCrossSectionIndustries.R	2013-08-09 11:32:24 UTC (rev 2749)
@@ -0,0 +1,52 @@
+
+##################################################################################################################
+#' This script computes the correlation between explicit, time-series industry factor returns and implicit, 
+#' cross-section industry factor returns, as 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_TimeSeriesVsCrossSectionIndustries.m"
+#'
+#' @author Xavier Valls \email{flamejat@@gmail.com}
+##################################################################################################################
+### Load data
+# loads weekly stock returns X and indices stock returns F
+load("../data/securitiesTS.Rda");
+Data_Securities = securitiesTS$data[ , -1 ]; # 1st column is date
+
+load("../data/sectorsTS.Rda");
+Data_Sectors = sectorsTS$data[ , -(1:2) ];
+
+load("../data/securitiesIndustryClassification.Rda");
+Securities_IndustryClassification = securitiesIndustryClassification$data;
+
+##################################################################################################################
+# Linear returns for stocks
+X = ( Data_Securities[ -1, ] - Data_Securities[ -nrow(Data_Securities),  ] ) / Data_Securities[ -nrow(Data_Securities),  ];
+
+# explicit, time-series industry factor returns
+F_ts = (Data_Sectors[ -1, ] - Data_Sectors[ -nrow(Data_Sectors),  ] ) / Data_Sectors[ -nrow(Data_Sectors),  ] ; 
+K = ncol(F_ts);
+
+# implicit, cross-section industry factor returns
+Sigma_X = (dim(X)[1]-1)/dim(X)[1] * cov(X);
+B    = Securities_IndustryClassification;
+Phi = diag(1 / diag( Sigma_X ), length(diag( Sigma_X ) ) );
+tmp = t(B) %*% Phi %*% B ;
+F_cs   = t( diag( diag(tmp) ^( -1 ), dim(tmp) ) %*% t(B) %*% Phi %*% t(X));
+
+
+##################################################################################################################
+### Correlation analysis
+Corr_cs_ts = matrix( 0, K, 1 );
+for( k in 1 : K )
+{
+    C = cor( cbind( F_cs[ , k], F_ts[ , k ] ));
+    Corr_cs_ts[ k ] = C[ 1, 2 ];
+}
+
+# time series factors are highly correlated with their cross-sectional counterparts
+dev.new();
+hist( Corr_cs_ts, seq( 0, 1, 0.1), xlim = c(-0.2, 0.2));
+