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

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

Author: xavierv
Date: 2013-08-05 21:22:22 +0200 (Mon, 05 Aug 2013)
New Revision: 2723

Added:
   pkg/Meucci/demo/S_FactorAnalysisNotOk.R
Log:
 - added FactorAnalysisNotOK demo script from chapter 3

Added: pkg/Meucci/demo/S_FactorAnalysisNotOk.R
===================================================================
--- pkg/Meucci/demo/S_FactorAnalysisNotOk.R	                        (rev 0)
+++ pkg/Meucci/demo/S_FactorAnalysisNotOk.R	2013-08-05 19:22:22 UTC (rev 2723)
@@ -0,0 +1,47 @@
+#'This script illustrates the hidden factor analysis puzzle, 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_FactorAnalysisNotOk.m"
+#'
+#' @author Xavier Valls \email{flamejat@@gmail.com}
+
+
+##################################################################################################################
+### Inputs
+
+N = 5; # market dimension
+K = 2; # factors dimension
+J = 10000; # numbers of simulations
+
+##################################################################################################################
+### Define true hidden loadings B 
+
+B = matrix( runif( N*K ), N, K ) - 0.5; 
+
+B = B / sqrt( 1.5 * max( max( B %*% t(B) ) ) );
+
+# define true hidden idiosyncratic variances 
+D = array( 1, N) - diag( B %*% t(B) );
+
+# define true hidden global covariance
+S = B %*% t( B ) + diag( D, length(D) );
+
+# generate normal variables with matching moments
+X = MvnRnd( matrix( 0, N, 1 ), S, J );
+
+# recover loadings FA$loadings, idiosyncratic variances FA$uniquenesness and factors FA$scores by factor analysis
+#[FA$loadings, FA$uniquenesness, T_, stats, F_] 
+FA = factanal(X, K, scores = "Bartlett" );
+
+# factor analysis recovers the structure exactly however...
+S_ = FA$loadings %*% t( FA$loadings ) + diag( FA$uniquenesses, length( FA$uniquenesses) );
+Match = 1 - max( abs( ( S - S_) / S) );
+print(Match);
+
+# ...the systematic+idiosyncratic decomposition is NOT recovered
+U_ = X - FA$scores %*% t(FA$loadings); # compute residuals
+S_U = cor( U_ ); # compute correlations
+# residuals are not idiosyncratic
+print( S_U );
\ No newline at end of file