[Returnanalytics-commits] r2851 - in pkg/Meucci: . R demo man

noreply at r-forge.r-project.org noreply at r-forge.r-project.org
Thu Aug 22 13:23:35 CEST 2013


Author: xavierv
Date: 2013-08-22 13:23:35 +0200 (Thu, 22 Aug 2013)
New Revision: 2851

Added:
   pkg/Meucci/R/MaxRsqCS.R
   pkg/Meucci/demo/S_CrossSectionConstrainedIndustries.R
   pkg/Meucci/man/MaxRsqCS.Rd
Modified:
   pkg/Meucci/DESCRIPTION
   pkg/Meucci/NAMESPACE
Log:
- added S_CrossSectionConstrainedIndustries demo script from chapter 3 and its associated functions

Modified: pkg/Meucci/DESCRIPTION
===================================================================
--- pkg/Meucci/DESCRIPTION	2013-08-22 10:48:35 UTC (rev 2850)
+++ pkg/Meucci/DESCRIPTION	2013-08-22 11:23:35 UTC (rev 2851)
@@ -90,3 +90,4 @@
     'CovertCompoundedReturns2Price.R'
     '
     FitOrnsteinUhlenbeck.R'
+    'MaxRsqCS.R'

Modified: pkg/Meucci/NAMESPACE
===================================================================
--- pkg/Meucci/NAMESPACE	2013-08-22 10:48:35 UTC (rev 2850)
+++ pkg/Meucci/NAMESPACE	2013-08-22 11:23:35 UTC (rev 2851)
@@ -25,6 +25,7 @@
 export(LognormalCopulaPdf)
 export(LognormalMoments2Parameters)
 export(LognormalParam2Statistics)
+export(MaxRsqCS)
 export(MleRecursionForStudentT)
 export(MvnRnd)
 export(NoisyObservations)

Added: pkg/Meucci/R/MaxRsqCS.R
===================================================================
--- pkg/Meucci/R/MaxRsqCS.R	                        (rev 0)
+++ pkg/Meucci/R/MaxRsqCS.R	2013-08-22 11:23:35 UTC (rev 2851)
@@ -0,0 +1,117 @@
+#' Solve for G that maximises sample r-square of X*G'*B' with X under constraints A*G<=D
+#' and Aeq*G=Deq (A,D, Aeq,Deq conformable matrices),as described in  A. Meucci, 
+#' "Risk and Asset Allocation", Springer, 2005.
+#'  
+#'  @param   X   : [matrix] (T x N)
+#'  @param   B   : [matrix] (T x K)
+#'  @param   W   : [matrix] (N x N)
+#'  @param   A   : [matrix] linear inequality constraints
+#'  @param   D   : [matrix]
+#'  @param   Aeq : [matrix] linear equality constraints
+#'  @param   Deq : [matrix] 
+#'  @param   lb  : [vector] lower bound
+#'  @param   ub  : [vector] upper bound
+#'  
+#'  @return   G   : [matrix] (N x K)
+#'
+#'  @note
+#'  Initial code by Tai-Ho Wang 
+#'
+#' @references
+#' \url{http://symmys.com/node/170}
+#' See Meucci's script for "MaxRsqCS.m"
+#'
+#' @author Xavier Valls \email{flamejat@@gmail.com}
+#' @export
+
+MaxRsqCS = function(X, B, W, A = NULL, D = NULL, Aeq = NULL, Deq, lb = NULL, ub = NULL)
+{
+    N = ncol(X);
+    K = ncol(B);
+
+    # compute sample estimates
+    Sigma_X = (dim(X)[1]-1)/dim(X)[1] * cov(X);
+
+
+    # restructure for feeding to quadprog 
+    Phi = t(W) %*% W;
+
+    # restructure the linear term of the objective function 
+    FirstDegree = matrix( Sigma_X %*% Phi %*% B, K * N, );
+
+    # restructure the quadratic term of the objective function
+    SecondDegree = Sigma_X;
+    
+    for( k in 1 : (N - 1) )
+    {
+        SecondDegree = blkdiag(SecondDegree, Sigma_X);
+    }
+    
+    SecondDegree = t( kron( sqrt(Phi) %*% B, diag( 1, N ) ) ) %*% SecondDegree %*% kron( sqrt( Phi ) %*% B, diag( 1, N ) );
+
+    # restructure the equality constraints 
+    if( !length(Aeq)  )
+    {
+        AEq = Aeq;
+    }else
+    {
+        AEq = blkdiag(Aeq);
+        for( k in 2 : K )
+        {
+            AEq = blkdiag(AEq, Aeq);
+        }
+    }
+
+    Deq = matrix( Deq, , 1);
+
+    # resturcture the inequality constraints 
+    if( length(A) )
+    {
+        AA = NULL
+        for( k in 1 : N )
+        {
+            AA = cbind( AA, kron(diag( 1, K ), A[ k ] ) ); ##ok<AGROW>
+        }
+    }else
+    {
+         AA = A;
+    }
+
+    if( length(D))
+    {
+        D = matrix( D, , 1 );
+    }
+
+    # restructure upper and lower bounds
+    if( length(lb) )
+    {
+        lb =  matrix( lb, K * N, 1 );
+    }
+
+    if( length(ub) )
+    {
+        ub = matrix( ub, K * N, 1 );
+    }
+
+    # initial guess
+    x0 = matrix( 1, K * N, 1 );
+    if(length(AA))
+    {
+        AA = ( AA + t(AA) ) / 2; # robustify
+          
+    }
+
+    Amat = rbind( AEq, AA);
+    bvec = c( Deq, D  );
+    
+    # solve the constrained generlized r-square problem by quadprog
+    #options = optimset('LargeScale', 'off', 'MaxIter', 2000, 'Display', 'none');
+    
+
+    b = ipop( c = matrix( FirstDegree ), H = SecondDegree, A = Amat, b = bvec, l = lb , u = ub , r = rep(0, length(bvec)) )
+    
+    # reshape for output
+    G = t( matrix( attributes(b)$primal, N, ) );
+
+    return( G );
+}
\ No newline at end of file

Added: pkg/Meucci/demo/S_CrossSectionConstrainedIndustries.R
===================================================================
--- pkg/Meucci/demo/S_CrossSectionConstrainedIndustries.R	                        (rev 0)
+++ pkg/Meucci/demo/S_CrossSectionConstrainedIndustries.R	2013-08-22 11:23:35 UTC (rev 2851)
@@ -0,0 +1,89 @@
+
+##################################################################################################################
+### This script fits a cross-sectional linear factor model creating industry factors, 
+### 
+### == Chapter 3 ==
+##################################################################################################################
+#' This script fits a cross-sectional linear factor model creating industry factors, where the industry factors 
+#' are constrained to be uncorrelated with the market 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_CrossSectionConstrainedIndustries.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/securitiesIndustryClassification.Rda");
+Securities_IndustryClassification = securitiesIndustryClassification$data;
+##################################################################################################################
+### Linear returns for stocks
+X = diff( Data_Securities ) / Data_Securities[ -nrow(Data_Securities), ];
+X = X[ ,1:30 ]; # consider first stocks only to lower run time (comment this line otherwise)
+
+T = dim(X)[1];
+N = dim(X)[2];
+B = Securities_IndustryClassification[ 1:N, ];
+K = dim(B)[ 2 ];
+m = array( 1/N, N );
+
+# compute sample estimates
+E_X = matrix( apply(X, 2, mean ) );
+Sigma_X = (dim(X)[1]-1)/dim(X)[1] * cov(X);
+
+# The optimal loadings turn out to be the standard multivariate weighted-OLS.
+Phi = diag(1 / diag( Sigma_X ), length(diag( Sigma_X ) ) );
+W 	= sqrt( Phi );
+
+##################################################################################################################
+### Solve for G that maximises sample r-square of X*G'*B' with X
+###  under constraints A*G<=D and Aeq*G=Deq (A,D, Aeq,Deq conformable matrices)
+A 	= NULL;
+Aeq = t(m) %*% t(Sigma_X);
+Deq = matrix( 0, K, 1 );
+#BOUNDARIES 
+lb = -100
+ub = 700
+G   = MaxRsqCS(X, B, W, A, D, Aeq, Deq, lb, ub);
+
+# compute intercept a and residual U
+F   = X %*% t(G);
+E_F = matrix(apply(F, 2, mean));
+a   = E_X - B %*% E_F;
+A_  = repmat(t(a), T, 1);
+U   = X - A_ - F %*% t(B);
+
+##################################################################################################################
+### Residual analysis
+
+M = cbind( U, F );
+SigmaJoint_UF = (dim(M)[1]-1)/dim(M)[1] * cov( M );
+CorrJoint_UF  = cov2cor(SigmaJoint_UF);
+Sigma_F = (dim(F)[1]-1)/dim(F)[1] * cov(F);
+Corr_F  = cov2cor( Sigma_F );
+Corr_F  = tril(Corr_F, -1);
+Corr_F  = Corr_F[ Corr_F != 0 ];
+
+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);
+
+Corr_UF = CorrJoint_UF[ 1:N, (N+1):(N+K) ];
+mean_Corr_UF = mean( abs(Corr_UF ) );
+max_Corr_UF  = max( abs(Corr_UF ) );
+disp(mean_Corr_U);
+disp(max_Corr_U);
+
+dev.new();
+hist(Corr_UF, 100);
\ No newline at end of file

Added: pkg/Meucci/man/MaxRsqCS.Rd
===================================================================
--- pkg/Meucci/man/MaxRsqCS.Rd	                        (rev 0)
+++ pkg/Meucci/man/MaxRsqCS.Rd	2013-08-22 11:23:35 UTC (rev 2851)
@@ -0,0 +1,48 @@
+\name{MaxRsqCS}
+\alias{MaxRsqCS}
+\title{Solve for G that maximises sample r-square of X*G'*B' with X under constraints A*G<=D
+and Aeq*G=Deq (A,D, Aeq,Deq conformable matrices),as described in  A. Meucci,
+"Risk and Asset Allocation", Springer, 2005.}
+\usage{
+  MaxRsqCS(X, B, W, A = NULL, D = NULL, Aeq = NULL, Deq,
+    lb = NULL, ub = NULL)
+}
+\arguments{
+  \item{X}{: [matrix] (T x N)}
+
+  \item{B}{: [matrix] (T x K)}
+
+  \item{W}{: [matrix] (N x N)}
+
+  \item{A}{: [matrix] linear inequality constraints}
+
+  \item{D}{: [matrix]}
+
+  \item{Aeq}{: [matrix] linear equality constraints}
+
+  \item{Deq}{: [matrix]}
+
+  \item{lb}{: [vector] lower bound}
+
+  \item{ub}{: [vector] upper bound}
+}
+\value{
+  G : [matrix] (N x K)
+}
+\description{
+  Solve for G that maximises sample r-square of X*G'*B'
+  with X under constraints A*G<=D and Aeq*G=Deq (A,D,
+  Aeq,Deq conformable matrices),as described in A. Meucci,
+  "Risk and Asset Allocation", Springer, 2005.
+}
+\note{
+  Initial code by Tai-Ho Wang
+}
+\author{
+  Xavier Valls \email{flamejat at gmail.com}
+}
+\references{
+  \url{http://symmys.com/node/170} See Meucci's script for
+  "MaxRsqCS.m"
+}
+



More information about the Returnanalytics-commits mailing list