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

noreply at r-forge.r-project.org noreply at r-forge.r-project.org
Wed Jul 31 18:08:57 CEST 2013


Author: xavierv
Date: 2013-07-31 18:08:56 +0200 (Wed, 31 Jul 2013)
New Revision: 2683

Added:
   pkg/Meucci/demo/S_ExactMeanAndCovariance.R
Modified:
   pkg/Meucci/DESCRIPTION
   pkg/Meucci/R/MvnRnd.R
   pkg/Meucci/man/MvnRnd.Rd
Log:
- added S_ExactMeanAndCovariance demo script and MvRnd function

Modified: pkg/Meucci/DESCRIPTION
===================================================================
--- pkg/Meucci/DESCRIPTION	2013-07-30 19:27:01 UTC (rev 2682)
+++ pkg/Meucci/DESCRIPTION	2013-07-31 16:08:56 UTC (rev 2683)
@@ -48,7 +48,8 @@
     latticeExtra,
     scatterplot3d,
     signal,
-    fExtremes
+    fExtremes,
+    QZ
 License: GPL
 URL: http://r-forge.r-project.org/projects/returnanalytics/
 Copyright: (c) 2012

Modified: pkg/Meucci/R/MvnRnd.R
===================================================================
--- pkg/Meucci/R/MvnRnd.R	2013-07-30 19:27:01 UTC (rev 2682)
+++ pkg/Meucci/R/MvnRnd.R	2013-07-31 16:08:56 UTC (rev 2683)
@@ -1,44 +1,47 @@
-#' Generates normal simulations whose sample moments match the population moments
+if ( !require( "QZ" ) ) stop("QZ package installation required for this script")
+
+#' Generate normal simulations whose sample moments match the population moments,
+#' as described in  A. Meucci, "Risk and Asset Allocation", Springer, 2005.
+#'  
+#'	@param   M : [vector] (N x 1) expectation
+#'	@param   S : [matrix] (N x N) covariance matrix
+#'	@param   J : [scalar] number of draws (even number)
+#'  
+#'	@return  X : [matrix] (J x N) of drawsF_U   : [vector] (J x 1) PDF values
 #'
-#' Adapted from file 'MvnRnd.m'. Most recent version of article and code available at http://www.symmys.com/node/162
-#' see A. Meucci - "Simulations with Exact Means and Covariances", Risk, July 2009
+#' @references
+#' \url{http://symmys.com/node/170}, \url{http://www.symmys.com/node/162}{A. Meucci - "Simulations with Exact Means and Covariances", Risk, July 2009}
+#' See Meucci's script for "MvnRnd.m"
 #'
-#' @param M         a numeric indicating the sample first moment of the distribution
-#' @param S         a covariance matrix
-#' @param J         a numeric indicating the number of trials
-#'
-#' @author Ram Ahluwalia \email{rahluwalia@@gmail.com}
-#' @references 
-#' \url{http://www.symmys.com}
-#' TODO: Add Schur decomposition. Right now function is only sampling from mvrnorm so sample moments do no match population moments
-#' I have sample code commented out below to implement this correctly but I require a function that returns the unitaryMatrix from a Schur decomposition
+#' @author Xavier Valls \email{flamejat@@gmail.com} and Ram Ahluwalia \email{rahluwalia@@gmail.com}
 #' @export
-MvnRnd = function( M , S , J )
+
+MvnRnd = function( M, S, J )
 {
-  library(MASS)
-  X = MASS::mvrnorm( n = J , mu = M , Sigma = S ) # Todo: need to swap with Meucci function and Schur method
-  return( X = X )
-    
-    # # compute sample covariance: NOTE defined as "cov(Y,1)", not as "cov(Y)"
-    # S_ = cov( Y )
-    # 
-    # # solve Riccati equation using Schur method
-    #     zerosMatrix = matrix( 0 ,  N ,  N );
-    #     # define the Hamiltonian matrix
-    #     H1 = cbind( zerosMatrix , -1*S_ )
-    #     H2 = cbind( -S , zerosMatrix ) 
-    #     H = rbind( H1 , H2 )
-    #     # perform its Schur decomposition. 
-    #     # TODO: check that the result returns real eigenvalues on the diagonal. ?Schur seems to give an example with real eigenvalues
-    #     schurDecomp = Schur( H )
-    #     T = SchurDecomp
-    #     # U_ = unitaryMatrix??? TODO: Find a function in R that returns the unitaryMatrix from a Schur decomposition
-    #     # U = ordschur(U_,T_,'lhp')
-    #     # U_lu = U(1:N,1:N)
-    #     # U_ld = U(N+1:end,1:N)
-    #     # B = U_ld/U_lu
-    # 
-    # # affine transformation to match mean and covariances
-    # # X = Y%*%B + repmat(M', J , 1 ) 
-    #
+	N = length(M);
+
+	# generate antithetic variables (mean = 0)
+	Y = rmvnorm( J/2, matrix( 0, N ), S );
+	Y = rbind( Y, -Y );
+
+	# compute sample covariance: NOTE defined as "cov(Y,1)", not as "cov(Y)"
+	S_ = ( dim(Y)[1] - 1 )/ dim(Y)[1] * cov( Y );
+
+	# solve Riccati equation using Schur method
+	H = rbind( cbind( matrix( 0, N, N ), -S ), cbind( -S, matrix( 0, N, N ) ) );
+	 
+	#Schur = Schur( H );
+	#U = ordschur(U_,T_,'lhp');
+	
+	U = ordqz( H, keyword = "lhp" )$Q;
+
+	U_lu = U[ 1:N, 1:N ];
+	U_ld = U[ (N+1):nrow(U), 1:N ];
+
+	B = U_ld %*% solve( U_lu );
+
+	# affine transformation to match mean and covariances
+	X = Y %*% B + kronecker( matrix( 1, J, 1 ),  t( M ) );
+
+	return( X );
 }
\ No newline at end of file

Added: pkg/Meucci/demo/S_ExactMeanAndCovariance.R
===================================================================
--- pkg/Meucci/demo/S_ExactMeanAndCovariance.R	                        (rev 0)
+++ pkg/Meucci/demo/S_ExactMeanAndCovariance.R	2013-07-31 16:08:56 UTC (rev 2683)
@@ -0,0 +1,41 @@
+#' Generate draws from a multivariate normal with matching mean and covariance, as described 
+#' in A. Meucci, "Risk and Asset Allocation", Springer, 2005,  Chapter 2.
+#'
+#' @references
+#' \url{http://symmys.com/node/170}
+#' See Meucci's script for "S_ExactMeanAndCovariance.m"
+#'
+#' @author Xavier Valls \email{flamejat@@gmail.com}
+
+
+########################################################################################################
+### Inputs
+N = 20;  # dimension (number of risk factors)
+J = 200; # number of simulations
+
+########################################################################################################
+### Generate desired population moments:
+
+# vector of expected values M
+M = matrix(runif( N ) -0.5);
+# covariance matrix S
+A  = matrix( runif( N * N ), c( N, N )) - 0.5;
+S = A %*% t( A );
+
+# generate sample of size J from multivariate normal N(M,S)
+#X = mvnrnd(M, S, J); # no match between sample and population moments (built-in) function
+X = MvnRnd( M, S, J ); # exact match between sample and population moments
+
+########################################################################################################
+### Compute sample moments and errors
+M_ = matrix( apply( X, 2, mean )); #apply
+S_ = ( dim( X )[1] - 1 )/ dim( X )[1] * cov( X );
+
+########################################################################################################
+### Check errors
+Err_M = max( abs( M - M_ ) ) / max( abs( M ) );
+Err_S = max( max( abs( S - S_) ) )/ max( max( abs( S ) ) );
+
+print(Err_M);
+print(Err_S);
+

Modified: pkg/Meucci/man/MvnRnd.Rd
===================================================================
--- pkg/Meucci/man/MvnRnd.Rd	2013-07-30 19:27:01 UTC (rev 2682)
+++ pkg/Meucci/man/MvnRnd.Rd	2013-07-31 16:08:56 UTC (rev 2683)
@@ -1,33 +1,34 @@
-\name{MvnRnd}
-\alias{MvnRnd}
-\title{Generates normal simulations whose sample moments match the population moments}
-\usage{
-  MvnRnd(M, S, J)
-}
-\arguments{
-  \item{M}{a numeric indicating the sample first moment of
-  the distribution}
-
-  \item{S}{a covariance matrix}
-
-  \item{J}{a numeric indicating the number of trials}
-}
-\description{
-  Adapted from file 'MvnRnd.m'. Most recent version of
-  article and code available at
-  http://www.symmys.com/node/162 see A. Meucci -
-  "Simulations with Exact Means and Covariances", Risk,
-  July 2009
-}
-\author{
-  Ram Ahluwalia \email{rahluwalia at gmail.com}
-}
-\references{
-  \url{http://www.symmys.com} TODO: Add Schur
-  decomposition. Right now function is only sampling from
-  mvrnorm so sample moments do no match population moments
-  I have sample code commented out below to implement this
-  correctly but I require a function that returns the
-  unitaryMatrix from a Schur decomposition
-}
-
+\name{MvnRnd}
+\alias{MvnRnd}
+\title{Generate normal simulations whose sample moments match the population moments,
+as described in  A. Meucci, "Risk and Asset Allocation", Springer, 2005.}
+\usage{
+  MvnRnd(M, S, J)
+}
+\arguments{
+  \item{M}{: [vector] (N x 1) expectation}
+
+  \item{S}{: [matrix] (N x N) covariance matrix}
+
+  \item{J}{: [scalar] number of draws (even number)}
+}
+\value{
+  X : [matrix] (J x N) of drawsF_U : [vector] (J x 1) PDF
+  values
+}
+\description{
+  Generate normal simulations whose sample moments match
+  the population moments, as described in A. Meucci, "Risk
+  and Asset Allocation", Springer, 2005.
+}
+\author{
+  Xavier Valls \email{flamejat at gmail.com} and Ram Ahluwalia
+  \email{rahluwalia at gmail.com}
+}
+\references{
+  \url{http://symmys.com/node/170},
+  \url{http://www.symmys.com/node/162}{A. Meucci -
+  "Simulations with Exact Means and Covariances", Risk,
+  July 2009} See Meucci's script for "MvnRnd.m"
+}
+



More information about the Returnanalytics-commits mailing list