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

noreply at r-forge.r-project.org noreply at r-forge.r-project.org
Fri Aug 30 10:41:18 CEST 2013


Author: xavierv
Date: 2013-08-30 10:41:17 +0200 (Fri, 30 Aug 2013)
New Revision: 2941

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

Modified: pkg/Meucci/DESCRIPTION
===================================================================
--- pkg/Meucci/DESCRIPTION	2013-08-30 08:39:38 UTC (rev 2940)
+++ pkg/Meucci/DESCRIPTION	2013-08-30 08:41:17 UTC (rev 2941)
@@ -96,7 +96,8 @@
     'FitOrnsteinUhlenbeck.R'
     'PlotVolVsCompositionEfficientFrontier.R'
     'BlackLittermanFormula.R'
+    'Log2Lin.R'
+    'PlotCompositionEfficientFrontier.R'
     '
     FitOrnsteinUhlenbeck.R'
-    'Log2Lin.R'
-    'PlotCompositionEfficientFrontier.R'
+    'MaxRsqTS.R'

Modified: pkg/Meucci/NAMESPACE
===================================================================
--- pkg/Meucci/NAMESPACE	2013-08-30 08:39:38 UTC (rev 2940)
+++ pkg/Meucci/NAMESPACE	2013-08-30 08:41:17 UTC (rev 2941)
@@ -31,6 +31,7 @@
 export(LognormalMoments2Parameters)
 export(LognormalParam2Statistics)
 export(MaxRsqCS)
+export(MaxRsqTS)
 export(MleRecursionForStudentT)
 export(MvnRnd)
 export(NoisyObservations)

Added: pkg/Meucci/R/MaxRsqTS.R
===================================================================
--- pkg/Meucci/R/MaxRsqTS.R	                        (rev 0)
+++ pkg/Meucci/R/MaxRsqTS.R	2013-08-30 08:41:17 UTC (rev 2941)
@@ -0,0 +1,126 @@
+#' Solve for B that maximises sample r-square of F'*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   F   : [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   B   : [matrix] (N x K)
+#'
+#'  @note
+#'  Initial code by Tai-Ho Wang 
+#'
+#' @references
+#' \url{http://symmys.com/node/170}
+#' See Meucci's script for "MaxRsqTS.m"
+#'
+#' @author Xavier Valls \email{flamejat@@gmail.com}
+#' @export
+
+MaxRsqTS = function(X, F, W, A = NULL, D = NULL, Aeq = NULL, Deq, lb = NULL, ub = NULL)
+{
+
+	N = dim(X)[ 2 ];
+	K = dim(F)[ 2 ];
+    X = matrix(as.numeric(X), dim(X))
+
+
+	# compute sample estimates
+	# E_X = apply( X, 2, mean);
+	# E_F = 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) ];
+
+
+    # restructure for feeding to quadprog 
+    Phi = t(W) %*% W;
+    trSigma_WX = sum( diag( Sigma_X %*% Phi ) );
+
+    # restructure the linear term of the objective function 
+    FirstDegree = ( -1 / trSigma_WX ) * matrix( t( Phi %*% Sigma_XF ), N*K, );
+
+    # restructure the quadratic term of the objective function
+    SecondDegree = Sigma_F;
+    for( k in 1 : (N - 1) )
+    {
+        SecondDegree = blkdiag( SecondDegree, Sigma_F );
+    }
+    SecondDegree = ( 1 / trSigma_WX ) * t(kron( sqrt( Phi ), diag( 1, K ))) %*% SecondDegree %*% kron( sqrt(Phi), diag( 1, K ) );
+
+    # restructure the equality constraints 
+    if( !length(Aeq)  )
+    {
+        AEq = Aeq;
+    }else
+    {
+        AEq = NULL;
+        for( k in 1 : N )
+        {
+            AEq = cbind( AEq, kron( diag( 1, K ), Aeq[k] ) );
+        }
+    }
+
+    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 = matrix( attributes(b)$primal, N, ) ;
+
+    return( G );
+}
\ No newline at end of file

Added: pkg/Meucci/demo/S_TimeSeriesConstrainedIndustries.R
===================================================================
--- pkg/Meucci/demo/S_TimeSeriesConstrainedIndustries.R	                        (rev 0)
+++ pkg/Meucci/demo/S_TimeSeriesConstrainedIndustries.R	2013-08-30 08:41:17 UTC (rev 2941)
@@ -0,0 +1,81 @@
+#' This script fits a time-series linear factor computing the industry factors loadings,  where the loadings are 
+#' bounded and constrained to yield unit exposure, 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_TimeSeriesConstrainedIndustries.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 date, 2nd column is SPX
+
+##################################################################################################################
+### Estimation
+# linear returns for stocks
+X = diff( Data_Securities ) / Data_Securities[ -nrow(Data_Securities), ];
+X = X[ ,1:20 ]; # consider first stocks only to lower run time (comment this line otherwise)
+
+# 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 ];
+
+##################################################################################################################
+# Solve for B that maximises sample r-square of F'*B' with X
+# under constraints A*B<=D and Aeq*B=Deq (A,D, Aeq,Deq conformable matrices)
+W   = diag( 1, N );
+A 	= NULL;
+D   = NULL;
+Aeq = matrix( 1, K, N ) / N
+Deq = matrix( 1, K, 1 );
+lb  = 0.8
+ub  = 1.2
+B   = MaxRsqTS(X, F, W, A, D, Aeq, Deq, lb, ub);
+
+# 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) ];
+
+# compute intercept a and residual U
+a = E_X - B %*% E_F;
+U = X - repmat(t(a), T, 1) - 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);
+
+# correlation 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));
+
+dev.new();
+hist(Corr_U, 100);
+
+# correlation between residuals and factors is not null
+Corr_UF = CorrJoint_UF[ 1:N, (N+1):(N+K) ];
+mean_Corr_UF = mean( abs(Corr_UF ) );
+max_Corr_UF  = max( abs(Corr_UF ) );
+
+dev.new();
+hist(Corr_UF, 100);
\ No newline at end of file

Added: pkg/Meucci/man/MaxRsqTS.Rd
===================================================================
--- pkg/Meucci/man/MaxRsqTS.Rd	                        (rev 0)
+++ pkg/Meucci/man/MaxRsqTS.Rd	2013-08-30 08:41:17 UTC (rev 2941)
@@ -0,0 +1,48 @@
+\name{MaxRsqTS}
+\alias{MaxRsqTS}
+\title{Solve for B that maximises sample r-square of F'*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{
+  MaxRsqTS(X, F, W, A = NULL, D = NULL, Aeq = NULL, Deq,
+    lb = NULL, ub = NULL)
+}
+\arguments{
+  \item{X}{: [matrix] (T x N)}
+
+  \item{F}{: [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{
+  B : [matrix] (N x K)
+}
+\description{
+  Solve for B that maximises sample r-square of F'*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
+  "MaxRsqTS.m"
+}
+



More information about the Returnanalytics-commits mailing list