[Returnanalytics-commits] r3953 - in pkg/Meucci: . R demo man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu Aug 13 18:11:52 CEST 2015
Author: xavierv
Date: 2015-08-13 18:11:51 +0200 (Thu, 13 Aug 2015)
New Revision: 3953
Modified:
pkg/Meucci/NAMESPACE
pkg/Meucci/R/MaxRsqCS.R
pkg/Meucci/R/MaxRsqTS.R
pkg/Meucci/R/MultivariateOUnCointegration.R
pkg/Meucci/demo/S_CPPI.R
pkg/Meucci/demo/S_CallsProjectionPricing.R
pkg/Meucci/demo/S_CheckDiagonalization.R
pkg/Meucci/demo/S_CornishFisher.R
pkg/Meucci/demo/S_CorrelationPriorUniform.R
pkg/Meucci/demo/S_CovarianceEvolution.R
pkg/Meucci/demo/S_CrossSectionConstrainedIndustries.R
pkg/Meucci/demo/S_CrossSectionIndustries.R
pkg/Meucci/demo/S_EigenvalueDispersion.R
pkg/Meucci/man/FitOU.Rd
pkg/Meucci/man/MaxRsqCS.Rd
pkg/Meucci/man/MaxRsqTS.Rd
pkg/Meucci/man/OUstep.Rd
pkg/Meucci/man/OUstepEuler.Rd
Log:
Formatted demo scripts and relating functions up to S_D*
Modified: pkg/Meucci/NAMESPACE
===================================================================
--- pkg/Meucci/NAMESPACE 2015-08-13 10:02:27 UTC (rev 3952)
+++ pkg/Meucci/NAMESPACE 2015-08-13 16:11:51 UTC (rev 3953)
@@ -27,6 +27,7 @@
export(Fit2Moms)
export(FitExpectationMaximization)
export(FitMultivariateGarch)
+export(FitOU)
export(FitOrnsteinUhlenbeck)
export(GenerateLogNormalDistribution)
export(GenerateUniformDrawsOnUnitSphere)
@@ -47,6 +48,8 @@
export(MvnRnd)
export(NoisyObservations)
export(NormalCopulaPdf)
+export(OUstep)
+export(OUstepEuler)
export(PHist)
export(PanicCopula)
export(PartialConfidencePosterior)
Modified: pkg/Meucci/R/MaxRsqCS.R
===================================================================
--- pkg/Meucci/R/MaxRsqCS.R 2015-08-13 10:02:27 UTC (rev 3952)
+++ pkg/Meucci/R/MaxRsqCS.R 2015-08-13 16:11:51 UTC (rev 3953)
@@ -1,9 +1,9 @@
-#' @title Solve for G that maximises sample r-square of X*G'*B' with X under constraints A*G<=D
-#' and Aeq*G=Deq
+#' @title Solve for G that maximises sample r-square of X*G'*B' with X under
+#' constraints A*G<=D and Aeq*G=Deq
#'
-#' @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.
+#' @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.
#'
#' @param X : [matrix] (T x N)
#' @param B : [matrix] (T x K)
@@ -17,105 +17,94 @@
#'
#' @return G : [matrix] (N x K)
#'
-#' @note
-#' Initial code by Tai-Ho Wang
+#' @note
+#' Initial code by Tai-Ho Wang
#'
#' @references
-#' A. Meucci - "Exercises in Advanced Risk and Portfolio Management" \url{http://symmys.com/node/170}.
-#' Used in "E 123 - Cross-section factors: generalized cross-section industry factors".
-#'
+#' A. Meucci - "Exercises in Advanced Risk and Portfolio Management"
+#' \url{http://symmys.com/node/170}. Used in "E 123 - Cross-section factors:
+#' generalized cross-section industry factors".
#' See Meucci's script for "MaxRsqCS.m"
#'
#' @author Xavier Valls \email{xaviervallspla@@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);
+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);
+ Sigma_X <- (dim(X)[1] - 1) / dim(X)[1] * cov(X)
- # restructure for feeding to quadprog
- Phi = t(W) %*% W;
+ # 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 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 <- 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);
+ 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);
+ 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>
+ # resturcture the inequality constraints
+ if(length(A)) {
+ AA <- NULL
+ for (k in 1 : N) {
+ AA <- cbind(AA, kron(diag(1, K), A[k]))
}
- }else
- {
- AA = A;
+ } else {
+ AA <- A
}
- if( length(D))
- {
- D = matrix( D, , 1 );
+ if (length(D)) {
+ D <- matrix(D, , 1)
}
# restructure upper and lower bounds
- if( length(lb) )
- {
- lb = matrix( lb, K * N, 1 );
+ if (length(lb)) {
+ lb <- matrix(lb, K * N, 1)
}
- if( length(ub) )
- {
- ub = matrix( ub, 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
-
+ x0 <- matrix(1, K * N, 1)
+ if(length(AA)) {
+ AA <- (AA + t(AA)) / 2 # robustify
}
- Amat = rbind( AEq, AA);
- bvec = c( Deq, D );
-
+ 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');
-
+ #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)) )
+ 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, ) );
+ G <- t(matrix(attributes(b)$primal, N, ))
- return( G );
-}
\ No newline at end of file
+ return(G)
+}
Modified: pkg/Meucci/R/MaxRsqTS.R
===================================================================
--- pkg/Meucci/R/MaxRsqTS.R 2015-08-13 10:02:27 UTC (rev 3952)
+++ pkg/Meucci/R/MaxRsqTS.R 2015-08-13 16:11:51 UTC (rev 3953)
@@ -1,6 +1,6 @@
-#' 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.
+#' 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)
@@ -18,109 +18,98 @@
#' Initial MATLAB's code by Tai-Ho Wang.
#'
#' @references
-#' A. Meucci - "Exercises in Advanced Risk and Portfolio Management" \url{http://symmys.com/node/170}.
-#' See Meucci's script for "MaxRsqTS.m"
+#' A. Meucci - "Exercises in Advanced Risk and Portfolio Management"
+#' \url{http://symmys.com/node/170}. See Meucci's script for "MaxRsqTS.m"
#'
#' @author Xavier Valls \email{xaviervallspla@@gmail.com}
#' @export
-MaxRsqTS = function(X, F, W, A = NULL, D = NULL, Aeq = NULL, Deq, lb = NULL, ub = NULL)
-{
+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))
+ 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) ];
+ # 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 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 <- 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 ) );
+ 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] ) );
+ # 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);
+ 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>
+ # resturcture the inequality constraints
+ if(length(A)) {
+ AA <- NULL
+ for (k in 1 : N) {
+ AA <- cbind(AA, kron(diag(1, K), A[k]))
}
- }else
- {
- AA = A;
+ } else {
+ AA <- A
}
- if( length(D))
- {
- D = matrix( D, , 1 );
+ if(length(D)) {
+ D <- matrix(D, , 1)
}
# restructure upper and lower bounds
- if( length(lb) )
- {
- lb = matrix( lb, K * N, 1 );
+ if(length(lb)) {
+ lb <- matrix(lb, K * N, 1)
}
- if( length(ub) )
- {
- ub = matrix( ub, 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
-
+ x0 <- matrix(1, K * N, 1)
+ if(length(AA)) {
+ AA <- (AA + t(AA)) / 2 # robustify
}
- Amat = rbind( AEq, AA);
- bvec = c( Deq, D );
+ 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');
-
+ #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)) )
-
+ 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, ) ;
+ G <- matrix(attributes(b)$primal, N,)
- return( G );
-}
\ No newline at end of file
+ return(G)
+}
Modified: pkg/Meucci/R/MultivariateOUnCointegration.R
===================================================================
--- pkg/Meucci/R/MultivariateOUnCointegration.R 2015-08-13 10:02:27 UTC (rev 3952)
+++ pkg/Meucci/R/MultivariateOUnCointegration.R 2015-08-13 16:11:51 UTC (rev 3953)
@@ -3,139 +3,158 @@
#' @param X_0 a matrix containing the starting value of each process
#' @param t a numeric containing the timestep
#' @param Mu a vector containing the unconditional expectation of the process
-#' @param Th a transition matrix, i.e., a fully generic square matrix that steers the deterministic portion
-#' of the evolution of the process
+#' @param Th a transition matrix, i.e., a fully generic square matrix that
+#' steers the deterministic portion of the evolution of the process
#' @param Sig a square matrix that drives the dispersion of the process
#'
#' @return a list containing
-#' @return X_t a vector containing the value of the process after the given timestep
+#' @return X_t a vector containing the value of the process after the given
+#' timestep
#' @return Mu_t a vector containing the conditional expectation of the process
#' @return Sig_t a matrix containing the covariance after the time step
#'
-#' \deqn{ X_{t+ \tau } = \big(I- e^{- \theta \tau } \big) \mu + e^{- \theta \tau } X_{t} + \epsilon _{t, \tau } }
+#' \deqn{ X_{t+ \tau } = \big(I- e^{- \theta \tau } \big) \mu +
+#' e^{- \theta \tau } X_{t} + \epsilon _{t, \tau } }
#' @references
-#' A. Meucci - "Review of Statistical Arbitrage, Cointegration, and Multivariate Ornstein-Uhlenbeck" - Formula (2)
-#' \url{http://ssrn.com/abstract=1404905}
+#' A. Meucci - "Review of Statistical Arbitrage, Cointegration, and Multivariate
+#' Ornstein-Uhlenbeck" - Formula (2) \url{http://ssrn.com/abstract=1404905}
+#'
#' @author Manan Shah \email{mkshah@@cmu.edu}
-OUstep = function( X_0 , t , Mu , Th , Sig )
-{
- NumSimul = nrow( X_0 )
- N = ncol( X_0 )
-
+#' @export
+
+OUstep <- function(X_0, t, Mu, Th, Sig) {
+ NumSimul <- nrow(X_0)
+ N <- ncol(X_0)
+
# location
- ExpM = expm( -Th * t )
-
+ ExpM <- expm(-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 ) + X_0 %*% ExpM
-
+ 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) +
+ 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( -TsT * t ) ) %*% VecSig
- Sig_t = VecSig_t
- dim( Sig_t ) = c( N , N )
- Sig_t = ( Sig_t + t( Sig_t ) ) / 2
-
- Eps = mvrnorm( NumSimul, rep( 0 , N ), Sig_t )
-
- X_t = Mu_t + Eps
- Mu_t = t( colMeans( Mu_t ) )
-
- return( list( X_t = X_t, Mu_t = Mu_t, Sig_t = Sig_t ) )
+ 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(-TsT * t)) %*% VecSig
+ Sig_t <- VecSig_t
+ dim(Sig_t) <- c(N, N)
+ Sig_t <- (Sig_t + t(Sig_t)) / 2
+
+ Eps <- mvrnorm(NumSimul, rep(0, N), Sig_t)
+
+ X_t <- Mu_t + Eps
+ Mu_t <- t(colMeans(Mu_t))
+
+ return(list(X_t = X_t, Mu_t = Mu_t, Sig_t = Sig_t))
}
-#' Generate the next element based on Ornstein-Uhlenbeck process using antithetic concept and assuming that the
-#' Brownian motion has Euler discretization
+#' Generate the next element based on Ornstein-Uhlenbeck process using
+#' antithetic concept and assuming that the Brownian motion has Euler
+#' discretization
#'
#' @param X_0 a matrix containing the starting value of each process
#' @param Dt a numeric containing the timestep
#' @param Mu a vector containing the unconditional expectation of the process
-#' @param Th a transition matrix, i.e., a fully generic square matrix that steers the deterministic portion
-#' of the evolution of the process
+#' @param Th a transition matrix, i.e., a fully generic square matrix that
+#' steers the deterministic portion of the evolution of the process
#' @param Sig a square matrix that drives the dispersion of the process
#'
#' @return a list containing
-#' @return X_t a vector containing the value of the process after the given timestep
+#' @return X_t a vector containing the value of the process after the given
+#' timestep
#' @return Mu_t a vector containing the conditional expectation of the process
#' @return Sig_t a matrix containing the covariance after the time step
#'
-#' \deqn{ X_{t+ \tau } = \big(I- e^{- \theta \tau } \big) \mu + e^{- \theta \tau } X_{t} + \epsilon _{t, \tau } }
+#' \deqn{ X_{t+ \tau } = \big(I- e^{- \theta \tau } \big) \mu +
+#' e^{- \theta \tau } X_{t} + \epsilon _{t, \tau } }
#' @references
-#' A. Meucci - "Review of Statistical Arbitrage, Cointegration, and Multivariate Ornstein-Uhlenbeck" - Formula (2)
-#' \url{http://ssrn.com/abstract=1404905}
+#' A. Meucci - "Review of Statistical Arbitrage, Cointegration, and Multivariate
+#' Ornstein-Uhlenbeck" - Formula (2). \url{http://ssrn.com/abstract=1404905}
+#'
#' @author Manan Shah \email{mkshah@@cmu.edu}
-OUstepEuler = function( X_0 , Dt , Mu , Th , Sig )
-{
- NumSimul = nrow( X_0 )
- N = ncol( X_0 )
-
+#' @export
+
+OUstepEuler <- function(X_0, Dt, Mu, Th, Sig){
+ NumSimul <- nrow(X_0)
+ N <- ncol(X_0)
+
# location
- ExpM = expm( as.matrix( -Th %*% Dt ) )
-
- # 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 ) + X_0 %*% ExpM
-
+ ExpM <- expm(as.matrix(-Th %*% Dt))
+
+ # 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)), mx * NumSimul, nx, byrow = T) +
+ X_0 %*% ExpM
+
# scatter
- Sig_t = Sig %*% Dt
- Eps = mvrnorm( NumSimul / 2, rep( 0 , N ) , Sig_t )
- Eps = rbind( Eps, -Eps)
-
- X_t = Mu_t + Eps
- Mu_t = t( colMeans( X_t ) )
-
- return( list( X_t = X_t, Mu_t = Mu_t, Sig_t = Sig_t ) )
+ Sig_t <- Sig %*% Dt
+ Eps <- mvrnorm(NumSimul / 2, rep(0, N), Sig_t)
+ Eps <- rbind(Eps, -Eps)
+
+ X_t <- Mu_t + Eps
+ Mu_t <- t(colMeans(X_t))
+
+ return(list(X_t = X_t, Mu_t = Mu_t, Sig_t = Sig_t))
}
-#' Fit the Ornstein-uhlenbeck process to model the behavior for different values of the timestep.
+#' Fit the Ornstein-uhlenbeck process to model the behavior for different values
+#' of the timestep.
#'
-#' @param Y a matrix containing the value of a process at various time steps.
+#' @param Y a matrix containing the value of a process at various time
+#' steps.
#' @param tau a numeric containing the timestep
#'
#' @return a list containing
#' @return Mu a vector containing the expectation of the process
-#' @return Sig a matrix containing the covariance of the resulting fitted OU process
-#' @return Th a transition matrix required for defining the fitted OU process
+#' @return Sig a matrix containing the covariance of the resulting fitted OU
+#' process
+#' @return Th a transition matrix required for defining the fitted OU
+#' process
#'
-#' \deqn{ x_{t+ \tau } = \big(I- e^{- \theta \tau } \big) \mu + e^{- \theta \tau } x_{t},
-#' vec \big( \Sigma _{ \tau } \big) \equiv \big( \Theta \oplus \Theta \big) ^{-1} \big(I- e^{( \Theta \oplus \Theta ) \tau } \big) vec \big( \Sigma \big) }
+#' \deqn{ x_{t+ \tau } = \big(I- e^{- \theta \tau } \big) \mu +
+#' e^{- \theta \tau } x_{t}, vec \big( \Sigma _{ \tau } \big) \equiv
+#' \big(\Theta \oplus \Theta \big) ^{-1} \big(I - e^{(\Theta \oplus \Theta)
+#' \tau } \big) vec \big( \Sigma \big) }
+#'
#' @references
-#' A. Meucci - "Review of Statistical Arbitrage, Cointegration, and Multivariate Ornstein-Uhlenbeck" - Formula (8),(9)
-#' \url{http://ssrn.com/abstract=1404905}
+#' A. Meucci - "Review of Statistical Arbitrage, Cointegration, and Multivariate
+#' Ornstein-Uhlenbeck" - Formula (8),(9). \url{http://ssrn.com/abstract=1404905}
#' @author Manan Shah \email{mkshah@@cmu.edu}
-FitOU = function ( Y, tau )
-{
+#' @export
+
+FitOU <- function(Y, tau) {
library(expm)
- T = nrow( Y )
- N = ncol( Y )
-
- X = Y[ -1 , ]
- F = cbind( rep( 1, T-1 ), Y [ 1:T-1 ,] )
- E_XF = t( X ) %*% F / T
- E_FF = t( F ) %*% F / T
- B = E_XF %*% solve( E_FF )
-
- Th = -logm ( B [ , -1 ] ) / tau
- Mu = solve( diag( N ) - B[ , -1 ] ) %*% B[ , 1 ]
-
- U = F %*% t( B ) - X
- Sig_tau = cov( U )
-
- N = length( Mu )
- TsT = kronecker( Th , diag( N ) ) + kronecker( diag( N ) , Th )
-
- VecSig_tau = Sig_tau
- dim( VecSig_tau ) = c( N^2 , 1 )
- VecSig = solve( diag( N^2 ) - expm( as.matrix( -TsT * tau ) ) ) %*% TsT %*% VecSig_tau
- Sig = VecSig
- dim( Sig ) = c( N , N )
-
- return( list( Mu = Mu, Th = Th, Sig = Sig ) )
-}
\ No newline at end of file
+ T <- nrow(Y)
+ N <- ncol(Y)
+
+ X <- Y[-1, ]
+ F <- cbind(rep(1, T - 1), Y [1:T - 1, ])
+ E_XF <- t(X) %*% F / T
+ E_FF <- t(F) %*% F / T
+ B <- E_XF %*% solve(E_FF)
+
+ Th <- -logm (B [, -1]) / tau
+ Mu <- solve(diag(N) - B[, -1]) %*% B[, 1]
+
+ U <- F %*% t(B) - X
+ Sig_tau <- cov(U)
+
+ N <- length(Mu)
+ TsT <- kronecker(Th, diag(N)) + kronecker(diag(N), Th)
+
+ VecSig_tau <- Sig_tau
+ dim(VecSig_tau) <- c(N ^ 2, 1)
+ VecSig <- solve(diag(N ^ 2) - expm(as.matrix(-TsT * tau))) %*% TsT %*%
+ VecSig_tau
+ Sig <- VecSig
+ dim(Sig) <- c(N, N)
+
+ return(list(Mu = Mu, Th = Th, Sig = Sig))
+}
Modified: pkg/Meucci/demo/S_CPPI.R
===================================================================
--- pkg/Meucci/demo/S_CPPI.R 2015-08-13 10:02:27 UTC (rev 3952)
+++ pkg/Meucci/demo/S_CPPI.R 2015-08-13 16:11:51 UTC (rev 3953)
@@ -1,114 +1,125 @@
-#' This script illustrates the CPPI (constant proportion portfolio insurance) dynamic strategy, as described in
-#' A. Meucci,"Risk and Asset Allocation", Springer, 2005, Chapter 6.
+#' This script illustrates the CPPI (constant proportion portfolio insurance)
+#' dynamic strategy, as described in A. Meucci,"Risk and Asset Allocation",
+#' Springer, 2005, Chapter 6.
#'
#' @references
-#' A. Meucci - "Exercises in Advanced Risk and Portfolio Management" \url{http://symmys.com/node/170},
-#' "E 264 - Constant proportion portfolio insurance".
+#' A. Meucci - "Exercises in Advanced Risk and Portfolio Management"
+#' \url{http://symmys.com/node/170}, "E 264 - Constant proportion portfolio
+#' insurance".
#'
#' See Meucci's script for "S_CPPI.m"
#
#' @author Xavier Valls \email{xaviervallspla@@gmail.com}
-##################################################################################################################
+################################################################################
### Input parameters
-Initial_Investment = 1000;
-Time_Horizon = 6 / 12; # in years
-Time_Step = 1 / 252; # in years
+Initial_Investment <- 1000
+Time_Horizon <- 6 / 12 # in years
+Time_Step <- 1 / 252 # in years
-m = 0.2; # yearly expected return on the underlying
-s = 0.40; # yearly expected percentage volatility on the stock index
-r = 0.04; # risk-free (money market) interest rate
+m <- 0.2 # yearly expected return on the underlying
+s <- 0.40 # yearly expected percentage volatility on the stock index
+r <- 0.04 # risk-free (money market) interest rate
-nSim = 30000;
+nSim <- 30000
-##################################################################################################################
+################################################################################
### Setup
# floor today (will evolve at the risk-free rate), e.g.: 950
-Floor = 980;
+Floor <- 980
# leverage on the cushion between your money and the floor, e.g. 3
-Multiple_CPPI = 5;
+Multiple_CPPI <- 5
-##################################################################################################################
+################################################################################
### Initialize values
-Underlying_Index = Initial_Investment; # value of the underlyting at starting time, normalzed to equal investment
-Start = Underlying_Index;
-Elapsed_Time = 0;
-Portfolio_Value = Initial_Investment;
+# value of the underlyting at starting time, normalzed to equal investment
+Underlying_Index <- Initial_Investment
+Start <- Underlying_Index
+Elapsed_Time <- 0
+Portfolio_Value <- Initial_Investment
-Cushion = max( 0, Portfolio_Value - Floor );
-Underlyings_in_Portfolio = min(Portfolio_Value, max( 0, Multiple_CPPI * Cushion ) );
-#Cash_in_Portfolio = Portfolio_Value - Underlyings_in_Portfolio;
-Underlying_in_Portfolio_Percent = Underlyings_in_Portfolio / Portfolio_Value;
+Cushion <- max(0, Portfolio_Value - Floor)
+Underlyings_in_Portfolio <- min(Portfolio_Value, max(0,
+ Multiple_CPPI * Cushion))
+#Cash_in_Portfolio <- Portfolio_Value - Underlyings_in_Portfolio
+Underlying_in_Portfolio_Percent <- Underlyings_in_Portfolio / Portfolio_Value
-Underlyings_in_Portfolio = Portfolio_Value * Underlying_in_Portfolio_Percent;
-Cash_in_Portfolio = Portfolio_Value - Underlyings_in_Portfolio;
+Underlyings_in_Portfolio <- Portfolio_Value * Underlying_in_Portfolio_Percent
+Cash_in_Portfolio <- Portfolio_Value - Underlyings_in_Portfolio
-##################################################################################################################
+################################################################################
### Initialize parameters for the plot (no theory in this)
-Portfolio_Series = Portfolio_Value;
-Market_Series = Underlying_Index;
-Percentage_Series = Underlying_in_Portfolio_Percent;
+Portfolio_Series <- Portfolio_Value
+Market_Series <- Underlying_Index
+Percentage_Series <- Underlying_in_Portfolio_Percent
-##################################################################################################################
+################################################################################
### Asset evolution and portfolio rebalancing
-while( Elapsed_Time < (Time_Horizon - 10^(-5)) ) # add this term to avoid errors
-{
+
+while (Elapsed_Time < (Time_Horizon - 10 ^ (-5)) ) {
# time elapses...
- Elapsed_Time = Elapsed_Time + Time_Step;
-
+ Elapsed_Time <- Elapsed_Time + Time_Step
+
# ...asset prices evolve and portfolio takes on new value...
- Multiplicator = exp( (m - s ^ 2 / 2) * Time_Step + s * sqrt( Time_Step ) * rnorm( NumSimul ));
- Underlying_Index = Underlying_Index * Multiplicator;
- Underlyings_in_Portfolio = Underlyings_in_Portfolio * Multiplicator;
- Cash_in_Portfolio = Cash_in_Portfolio * exp(r * Time_Step);
- Portfolio_Value = Underlyings_in_Portfolio + Cash_in_Portfolio;
+ Multiplicator <- exp((m - s ^ 2 / 2) * Time_Step + s *
+ sqrt(Time_Step) * rnorm(NumSimul))
+ Underlying_Index <- Underlying_Index * Multiplicator
+ Underlyings_in_Portfolio <- Underlyings_in_Portfolio * Multiplicator
+ Cash_in_Portfolio <- Cash_in_Portfolio * exp(r * Time_Step)
+ Portfolio_Value <- Underlyings_in_Portfolio + Cash_in_Portfolio
# ...and we rebalance our portfolio
- Floor = Floor * exp( r * Time_Step );
- Cushion = pmax( 0, (Portfolio_Value - Floor) );
- Underlyings_in_Portfolio = pmin(Portfolio_Value, pmax( 0, Multiple_CPPI * Cushion) );
- Cash_in_Portfolio = Portfolio_Value - Underlyings_in_Portfolio;
- Underlying_in_Portfolio_Percent = Underlyings_in_Portfolio / Portfolio_Value;
-
+ Floor <- Floor * exp(r * Time_Step)
+ Cushion <- pmax(0, (Portfolio_Value - Floor))
+ Underlyings_in_Portfolio <- pmin(Portfolio_Value, pmax(0,
+ Multiple_CPPI * Cushion))
+ Cash_in_Portfolio <- Portfolio_Value -
+ Underlyings_in_Portfolio
+ Underlying_in_Portfolio_Percent <- Underlyings_in_Portfolio /
+ Portfolio_Value
+
# store one path for the movie (no theory in this)
- Portfolio_Series = cbind( Portfolio_Series, Portfolio_Value[ 1 ] ); ##ok<*AGROW>
- Market_Series = cbind( Market_Series, Underlying_Index[ 1 ] );
- Percentage_Series = cbind( Percentage_Series, Underlying_in_Portfolio_Percent[ 1 ] );
+ Portfolio_Series <- cbind(Portfolio_Series, Portfolio_Value[1])
+ Market_Series <- cbind(Market_Series, Underlying_Index[1])
+ Percentage_Series <- cbind(Percentage_Series,
+ Underlying_in_Portfolio_Percent[1])
}
-##################################################################################################################
+################################################################################
### Play the movie for one path
-Time = seq( 0, Time_Horizon, Time_Step);
-y_max = max( cbind( Portfolio_Series, Market_Series) ) * 1.2;
-dev.new();
-par( mfrow = c(2,1))
-for( i in 1 : length(Time) )
-{
- plot( Time[ 1:i ], Portfolio_Series[ 1:i ], type ="l", lwd = 2.5, col = "blue", ylab = "value",
- xlim = c(0, Time_Horizon), ylim = c(0, y_max), main = "investment (blue) vs underlying (red) value");
- lines( Time[ 1:i ], Market_Series[ 1:i ], lwd = 2, col = "red" );
- #axis( 1, [0, Time_Horizon, 0, y_max]);
-
- plot(Time[ 1:i ], Percentage_Series[ 1:i ], type = "h", col = "red", xlab = "time", ylab = "#",
- xlim = c(0, Time_Horizon), ylim =c(0,1), main = "percentage of underlying in portfolio");
-
+Time <- seq(0, Time_Horizon, Time_Step)
+y_max <- max(cbind(Portfolio_Series, Market_Series)) * 1.2
+dev.new()
+par(mfrow <- c(2,1))
+for (i in 1 : length(Time)) {
+ plot(Time[1:i], Portfolio_Series[1:i], type ="l", lwd = 2.5,
+ col = "blue", ylab = "value", xlim = c(0, Time_Horizon),
+ ylim = c(0, y_max),
+ main = "investment (blue) vs underlying (red) value")
+ lines(Time[1:i], Market_Series[1:i], lwd = 2, col = "red")
+ #axis(1, [0, Time_Horizon, 0, y_max])
+
+ plot(Time[1:i], Percentage_Series[1:i], type = "h", col = "red",
+ xlab = "time", ylab = "#", xlim = c(0, Time_Horizon), ylim =c(0,1),
+ main = "percentage of underlying in portfolio")
}
-
-##################################################################################################################
+################################################################################
### plot the scatterplot
-dev.new();
+dev.new()
# marginals
-NumBins = round(10 * log(NumSimul));
-layout( matrix(c(1,2,2,2,1,2,2,2,1,2,2,2,0,3,3,3), 4, 4, byrow = TRUE));
-barplot( table( cut( Portfolio_Value, NumBins )), horiz=TRUE, yaxt="n");
+NumBins <- round(10 * log(NumSimul))
+layout(matrix(c(1, 2, 2, 2, 1, 2, 2, 2, 1, 2, 2, 2, 0, 3, 3, 3), 4, 4,
+ byrow = TRUE))
+barplot(table(cut(Portfolio_Value, NumBins)), horiz = TRUE, yaxt = "n")
# joint scatter plot
-plot(Underlying_Index, Portfolio_Value, xlab = "underlying at horizon (~ buy & hold )", ylab = "investment at horizon" );
-so = sort( Underlying_Index );
-lines( so, so, col = "red" );
+plot(Underlying_Index, Portfolio_Value,
+ xlab = "underlying at horizon (~ buy & hold)",
+ ylab = "investment at horizon")
+so <- sort(Underlying_Index)
+lines(so, so, col = "red")
-barplot( table( cut( Underlying_Index, NumBins )), yaxt="n");
-
+barplot(table(cut(Underlying_Index, NumBins)), yaxt = "n")
Modified: pkg/Meucci/demo/S_CallsProjectionPricing.R
===================================================================
--- pkg/Meucci/demo/S_CallsProjectionPricing.R 2015-08-13 10:02:27 UTC (rev 3952)
+++ pkg/Meucci/demo/S_CallsProjectionPricing.R 2015-08-13 16:11:51 UTC (rev 3953)
@@ -1,104 +1,111 @@
-#'This script projects the distribution of the market invariants for the derivatives market
-#'Then it computes the distribution of prices at the investment horizon as described in A. Meucci,
-#'"Risk and Asset Allocation", Springer, 2005, Chapter 3.
+#'This script projects the distribution of the market invariants for the
+#' derivatives market
+#'Then it computes the distribution of prices at the investment horizon as
+#' described in A. Meucci, "Risk and Asset Allocation", Springer, 2005, Ch 3.
#'
#' @references
-#' A. Meucci - "Exercises in Advanced Risk and Portfolio Management" \url{http://symmys.com/node/170},
-#' "E 143 - Derivatives market: projection of invariants".
+#' A. Meucci - "Exercises in Advanced Risk and Portfolio Management"
+#' \url{http://symmys.com/node/170}, "E 143 - Derivatives market: projection of
+#' invariants".
#'
#' See Meucci's script for "S_CallsProjectionPricing.m"
#'
#' @author Xavier Valls \email{xaviervallspla@@gmail.com}
-##################################################################################################################
+################################################################################
### Load data
# load 'spot' for underlying and current vol surface, given by
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/returnanalytics -r 3953
More information about the Returnanalytics-commits
mailing list