[Returnanalytics-commits] r2242 - in pkg/PerformanceAnalytics/sandbox/Meucci: R data demo
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sun Aug 19 02:36:14 CEST 2012
Author: mkshah
Date: 2012-08-19 02:36:14 +0200 (Sun, 19 Aug 2012)
New Revision: 2242
Added:
pkg/PerformanceAnalytics/sandbox/Meucci/data/DB_SwapParRates.rda
pkg/PerformanceAnalytics/sandbox/Meucci/demo/S_FitProjectRates.R
Modified:
pkg/PerformanceAnalytics/sandbox/Meucci/R/EmpiricalMultivariateOUnCointegration.R
pkg/PerformanceAnalytics/sandbox/Meucci/R/MeanDiversificationFrontier.R
pkg/PerformanceAnalytics/sandbox/Meucci/data/00index
pkg/PerformanceAnalytics/sandbox/Meucci/demo/MeanDiversificationFrontier.R
Log:
Correcting functions, adding demos and data files
Modified: pkg/PerformanceAnalytics/sandbox/Meucci/R/EmpiricalMultivariateOUnCointegration.R
===================================================================
--- pkg/PerformanceAnalytics/sandbox/Meucci/R/EmpiricalMultivariateOUnCointegration.R 2012-08-18 16:40:50 UTC (rev 2241)
+++ pkg/PerformanceAnalytics/sandbox/Meucci/R/EmpiricalMultivariateOUnCointegration.R 2012-08-19 00:36:14 UTC (rev 2242)
@@ -1,5 +1,6 @@
FitOU = function ( Y, tau )
{
+ library(expm)
T = nrow( Y )
N = ncol( Y )
@@ -20,7 +21,7 @@
VecSig_tau = Sig_tau
dim( VecSig_tau ) = c( N^2 , 1 )
- VecSig = solve( diag( N^2 ) - expm( Matrix( -TsT * tau ) ) ) %*% TsT %*% VecSig_tau
+ VecSig = solve( diag( N^2 ) - expm( as.matrix( -TsT * tau ) ) ) %*% TsT %*% VecSig_tau
Sig = VecSig
dim( Sig ) = c( N , N )
@@ -33,20 +34,20 @@
N = ncol( X_0 )
# location
- ExpM = expm( Matrix( -Th %*% t ) )
+ ExpM = expm( as.matrix ( -Th * t ) )
# repmat = function(X,m,n) - R equivalent of repmat (matlab)
X = t( Mu - ExpM %*% Mu )
mx = dim( X )[1]
nx = dim( X )[2]
- Mu_t = matrix( t ( matrix( X , mx , nx*1 ) ), mx * NumSimul, nx * 1, byrow = T ) + t( X_0 %*% ExpM )
+ Mu_t = matrix( t ( matrix( X , mx , nx*1 ) ), mx * NumSimul, nx * 1, byrow = T ) + X_0 %*% ExpM
# scatter
TsT = kronecker( Th , diag( N ) ) + kronecker( diag( N ) , Th )
VecSig = Sig
dim( VecSig ) = c( N^2 , 1 )
- VecSig_t = solve( TsT ) %*% ( diag( N^2 ) - expm( Matrix( -TsT %*% t ) ) ) %*% VecSig
+ VecSig_t = solve( TsT ) %*% ( diag( N^2 ) - expm( as.matrix( -TsT * t ) ) ) %*% VecSig
Sig_t = VecSig_t
dim( Sig_t ) = c( N , N )
Sig_t = ( Sig_t + t( Sig_t ) ) / 2
@@ -64,14 +65,14 @@
N = length( x_0 )
# location
- Mu_t = Mu + expm( Matrix( -Th %*% t ) ) %*% ( x_0 - Mu )
+ Mu_t = Mu + expm( as.matrix( -Th * t ) ) %*% ( x_0 - Mu )
# scatter
TsT = kronecker( Th , diag( N ) ) + kronecker( diag( N ) , Th )
VecSig = Sig
dim( VecSig ) = c( N^2 , 1 )
- VecSig_t = solve( TsT ) %*% ( diag( N^2 ) - expm( Matrix( -TsT %*% t ) ) ) %*% VecSig
+ VecSig_t = solve( TsT ) %*% ( diag( N^2 ) - expm( as.matrix( -TsT * t ) ) ) %*% VecSig
Sig_t = VecSig_t
dim( Sig_t ) = c( N , N )
}
\ No newline at end of file
Modified: pkg/PerformanceAnalytics/sandbox/Meucci/R/MeanDiversificationFrontier.R
===================================================================
--- pkg/PerformanceAnalytics/sandbox/Meucci/R/MeanDiversificationFrontier.R 2012-08-18 16:40:50 UTC (rev 2241)
+++ pkg/PerformanceAnalytics/sandbox/Meucci/R/MeanDiversificationFrontier.R 2012-08-19 00:36:14 UTC (rev 2242)
@@ -51,20 +51,20 @@
e = GenFirstEigVect( S , B )
E = cbind( E , e )
}
-
+
for ( n in N - K + 1:N )
{
B = t( E ) %*% S
e = GenFirstEigVect( S , B )
E = cbind( E , e )
}
-
+
# swap order
E = cbind( E[ , N - K + 1:N ], E[ , 1:N - K ] )
- }
-
- v = t( E ) %*% S %*% E
- L = diag( v , nrow = length( v ) )
+ }
+
+ v = t( E ) %*% as.matrix( S ) %*% E
+ L = diag( v )
G = diag( sqrt( L ), nrow = length( L ) ) %*% solve( E )
G = G[ K + 1:N , ]
@@ -73,12 +73,13 @@
MaxEntropy = function( G , w_b , w_0 , Constr )
{
+ library( nloptr )
# Nested function that computes fitness
nestedfun = function( x )
{
- v_ = G %*% ( x - w_b )
+ v_ = G %*% ( x - as.matrix( w_b ) )
p = v_ * v_
- R_2 = max( 10^(-10), p / colSums( p ) )
+ R_2 = pmax( 10^(-10), p / colSums( p ) )
Minus_Ent = t( R_2 ) * log( R_2 )
# evaluate gradient
@@ -88,12 +89,13 @@
}
local_opts <- list( algorithm = "NLOPT_LD_SLSQP", xtol_rel = 1.0e-6 ,
- check_derivatives = TRUE , check_derivatives_print = "all" ,
+ check_derivatives = FALSE , #check_derivatives_print = "all" ,
eval_f = nestedfun )
- x = nloptr( x0 = x0 , eval_f = nestedfunC ,
+ x = nloptr( x0 = w_0 , eval_f = nestedfun ,
opts = list( algorithm = "NLOPT_LD_AUGLAG" , local_opts = local_opts ,
print_level = 2 , maxeval = 1000 ,
- check_derivatives = TRUE , check_derivatives_print = "all" , xtol_rel = 1.0e-6 ) )
+ check_derivatives = FALSE , #check_derivatives_print = "all" ,
+ xtol_rel = 1.0e-6 ) )
return( x )
}
@@ -105,10 +107,10 @@
# compute frontier extrema]
library( limSolve )
- w_MaxExp = linp( E = Constr$Aeq , F = Constr$beq , G = -1*Constr$A , H = -1*Constr$b, Cost = -Mu , ispos = FALSE)$X
- MaxExp = t( Mu ) %*% ( w_MaxExp - w_b )
+ w_MaxExp = as.matrix( linp( E = Constr$Aeq , F = Constr$beq , G = -1*Constr$A , H = -1*Constr$b, Cost = -t( Mu ) , ispos = FALSE)$X )
+ MaxExp = t( Mu ) %*% ( w_MaxExp - as.matrix( w_b ) )
- w_MaxNe = MaxEntropy( G , w_b , w_0 , Constr )
+ w_MaxNe = MaxEntropy( GenPCBasisResult$G , w_b , w_0 , Constr )
ExpMaxNe = t( Mu ) %*% ( w_MaxNe - w_b )
# slice efficient frontier in NumPortf equally thick horizontal sections
@@ -131,13 +133,13 @@
ConstR$Aeq = cbind( Constr$Aeq, t( Mu ) )
ConstR$beq = cbind( Constr$beq, TargetExp[ k ] + t( Mu ) %*% w_b )
- w = MaxEntropy( G , w_b , w_0 , ConstR )
+ w = MaxEntropy( GenPCBasisResult$G , w_b , w_0 , ConstR )
m = t( Mu ) %*% ( w - w_b )
s = sqrt( t( w - w_b ) %*% S %*% ( w - w_b ) )
- v_ = G %*% ( w - w_b )
+ v_ = GenPCBasisResult$G %*% ( w - w_b )
TE_contr = v_ * v_ / s
R_2 = max( 10^(-10) , TE_contr / colSums( TE_contr ) )
Modified: pkg/PerformanceAnalytics/sandbox/Meucci/data/00index
===================================================================
--- pkg/PerformanceAnalytics/sandbox/Meucci/data/00index 2012-08-18 16:40:50 UTC (rev 2241)
+++ pkg/PerformanceAnalytics/sandbox/Meucci/data/00index 2012-08-19 00:36:14 UTC (rev 2242)
@@ -5,4 +5,6 @@
MeucciTweakTest A_ Aeq_ b_ beq_ db_ g_ lb_ ub_
pseudodata data
ReturnsDistribution P X
-SectorsSnP500 DP P
\ No newline at end of file
+SectorsSnP500 DP P
+MeanDiversificationFrontier S Mu w_b
+DB_SwapParRates Rates Dates
\ No newline at end of file
Added: pkg/PerformanceAnalytics/sandbox/Meucci/data/DB_SwapParRates.rda
===================================================================
(Binary files differ)
Property changes on: pkg/PerformanceAnalytics/sandbox/Meucci/data/DB_SwapParRates.rda
___________________________________________________________________
Added: svn:mime-type
+ application/octet-stream
Modified: pkg/PerformanceAnalytics/sandbox/Meucci/demo/MeanDiversificationFrontier.R
===================================================================
--- pkg/PerformanceAnalytics/sandbox/Meucci/demo/MeanDiversificationFrontier.R 2012-08-18 16:40:50 UTC (rev 2241)
+++ pkg/PerformanceAnalytics/sandbox/Meucci/demo/MeanDiversificationFrontier.R 2012-08-19 00:36:14 UTC (rev 2242)
@@ -16,9 +16,9 @@
# long-short constraints...
Constr = list()
Constr$A = rbind( diag( N ), -diag( N ) )
-Constr$b = rbind( rep( 1, N ), rep( 0.1, N ) )
-Constr$Aeq = rep( 1 , N ) # budget constraint...
-Constr$beq = 1
+Constr$b = rbind( as.matrix( rep( 1, N ) ), as.matrix( rep( 0.1, N ) ) )
+Constr$Aeq = t( as.matrix( rep( 1 , N ) ) ) # budget constraint...
+Constr$beq = as.matrix( 1 )
# mean-diversification analysis and frontier
EntropyFrontier = MeanTCEntropyFrontier( S , Mu , w_b , w_0 , Constr )
Added: pkg/PerformanceAnalytics/sandbox/Meucci/demo/S_FitProjectRates.R
===================================================================
--- pkg/PerformanceAnalytics/sandbox/Meucci/demo/S_FitProjectRates.R (rev 0)
+++ pkg/PerformanceAnalytics/sandbox/Meucci/demo/S_FitProjectRates.R 2012-08-19 00:36:14 UTC (rev 2242)
@@ -0,0 +1,38 @@
+# This script fits the swap rates dynamics to a multivariate Ornstein-Uhlenbeck process
+# and computes and plots the estimated future distribution
+# see A. Meucci (2009)
+# "Review of Statistical Arbitrage, Cointegration, and Multivariate Ornstein-Uhlenbeck"
+# available at ssrn.com
+
+# Code by A. Meucci, April 2009
+# Most recent version available at www.symmys.com > Teaching > MATLAB
+
+# inputs
+TimeStep = 5 # select time interval (days)
+Taus = c( 1/252,5/252,1/12,.5,1,2,10 ) # select horizon projection (years)
+Pick = c( 2,3 )
+
+# estimation
+StepRates = Rates[ seq( from = 1, to = nrow( Rates ), by = TimeStep ) , ]
+OUResult = FitOU( as.matrix( StepRates ) , TimeStep / 252 )
+
+for( s in 1:length(Taus))
+{
+ RGB = .6 * as.matrix( c( runif(1),runif(1),runif(1) ) )
+ tau = Taus[ s ]
+
+ # projection
+ # x_T = Mu
+ x_T = Rates[ nrow( Rates ) , ]
+ OUstepResult = OUstep( as.matrix( x_T ) , tau , OUResult$Mu , OUResult$Th , OUResult$Sig )
+
+ # plot
+ # historical observations
+ plot( Rates[ , Pick[1] ] , Rates[ , Pick[2] ] )
+
+ # current observation
+ plot( x_T[ Pick[1] ] , x_T[ Pick[2] ] )
+
+ # horizon location
+ plot( OUstepResult$Mu_t[ Pick[1] ] , OUstepResult$Mu_t[ Pick[2] ] )
+}
\ No newline at end of file
More information about the Returnanalytics-commits
mailing list