[Returnanalytics-commits] r2662 - in pkg/Meucci: . demo man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Mon Jul 29 14:18:32 CEST 2013
Author: xavierv
Date: 2013-07-29 14:18:31 +0200 (Mon, 29 Jul 2013)
New Revision: 2662
Added:
pkg/Meucci/demo/S_EstimateMomentsComboEvaluation.R
pkg/Meucci/demo/S_Toeplitz.R
pkg/Meucci/demo/S_VolatilityClustering.R
Modified:
pkg/Meucci/DESCRIPTION
pkg/Meucci/demo/S_EstimateExpectedValueEvaluation.R
pkg/Meucci/man/Central2Raw.Rd
pkg/Meucci/man/Cumul2Raw.Rd
pkg/Meucci/man/Raw2Central.Rd
pkg/Meucci/man/Raw2Cumul.Rd
Log:
-added two demo scripts from chapter 3 and one from chapter 4. Fixed documentation
Modified: pkg/Meucci/DESCRIPTION
===================================================================
--- pkg/Meucci/DESCRIPTION 2013-07-29 10:37:59 UTC (rev 2661)
+++ pkg/Meucci/DESCRIPTION 2013-07-29 12:18:31 UTC (rev 2662)
@@ -76,11 +76,7 @@
'SimulateJumpDiffusionMerton.R'
'BlackScholesCallPrice.R'
'InterExtrapolate.R'
- 'Central2Raw.R'
'CentralAndStandardizedStatistics.R'
- 'Cumul2Raw.R'
- 'Raw2Central.R'
- 'Raw2Cumul.R'
'FitExpectationMaximization.R'
'QuantileMixture.R'
'GenerateUniformDrawsOnUnitSphere.R'
Modified: pkg/Meucci/demo/S_EstimateExpectedValueEvaluation.R
===================================================================
--- pkg/Meucci/demo/S_EstimateExpectedValueEvaluation.R 2013-07-29 10:37:59 UTC (rev 2661)
+++ pkg/Meucci/demo/S_EstimateExpectedValueEvaluation.R 2013-07-29 12:18:31 UTC (rev 2662)
@@ -3,7 +3,7 @@
#'
#' @references
#' \url{http://symmys.com/node/170}
-#' See Meucci's script for "S_EigenValueprintersion.R"
+#' See Meucci's script for "S_EigenValueDispersion.R"
#'
#' @author Xavier Valls \email{flamejat@@gmail.com}
@@ -21,8 +21,8 @@
i_T = matrix( rlnorm( T, Mu, Sigma ), 1, T); # series generated by "nature": do not know the distribution
-G_Hat_1 = function(X) X[ , 1 ] * X[ ,3 ]; # estimator of unknown functional G_1=x(1)*x(3)
-G_Hat_2 = function(X) apply( X, 1,mean); # estimator of unknown functional G_1=sample mean
+G_Hat_1 = function(X) X[ , 1 ] * X[ , 3 ]; # estimator of unknown functional G_1=x(1)*x(3)
+G_Hat_2 = function(X) apply( X, 1, mean ); # estimator of unknown functional G_1=sample mean
G1 = G_Hat_1( i_T );
G2 = G_Hat_2( i_T );
@@ -59,19 +59,19 @@
par( mfrow = c(2,1) );
hist(G1, NumBins);
points(G_fX, 0, pch = 21, bg = "red", main = "estimator: x(1)*x(3)");
-#set(h, 'markersize', 20, 'col', 'r');
+
hist(G2, NumBins);
points(G_fX, 0, pch = 21, bg = "red", main = "estimator: sample mean" );
-#set(h, 'markersize', 20, 'col', 'r');
+
# loss
dev.new();
par( mfrow = c(2,1) );
-hist(Loss_G1, NumBins, main = "estimator: x(1)*x(3)");
+hist(Loss_G1, NumBins, main = "loss of estimator: x(1)*x(3)");
-hist(Loss_G2, NumBins, main = "estimator: sample mean" );
+hist(Loss_G2, NumBins, main = "loss of estimator: sample mean" );
##################################################################################################################
@@ -116,13 +116,13 @@
dev.new();
NumBins = round(10*log(nSim));
- par( mfrow = c(2,1) );
-
- hist(G1, NumBins);
- points(G_fX, 0, pch = 21, bg = "red", main = "estimator: x(1)*x(3)");
-
- hist(G2, NumBins);
- points(G_fX, 0, pch = 21, bg = "red", main = "estimator: sample mean" );
+ par( mfrow = c(2,1) );
+
+ hist(G1, NumBins);
+ points(G_fX, 0, pch = 21, bg = "red", main = "estimator: x(1)*x(3)");
+
+ hist(G2, NumBins);
+ points(G_fX, 0, pch = 21, bg = "red", main = "estimator: sample mean" );
}
Added: pkg/Meucci/demo/S_EstimateMomentsComboEvaluation.R
===================================================================
--- pkg/Meucci/demo/S_EstimateMomentsComboEvaluation.R (rev 0)
+++ pkg/Meucci/demo/S_EstimateMomentsComboEvaluation.R 2013-07-29 12:18:31 UTC (rev 2662)
@@ -0,0 +1,209 @@
+#'This script familiarizes the user with the evaluation of an estimator:replicability, loss, error,
+#'bias and inefficiency as described in A. Meucci,"Risk and Asset Allocation", Springer, 2005, Chapter 4.
+#'
+#' @references
+#' \url{http://symmys.com/node/170}
+#' See Meucci's script for "S_EstimateMomentsComboEvaluation.m"
+#'
+#' @author Xavier Valls \email{flamejat@@gmail.com}
+
+##################################################################################################################
+### Set parametesr
+
+T = 10;
+a = 0.5;
+m_Y = 0.1;
+s_Y = 0.2;
+m_Z = 0;
+s_Z = 0.15;
+
+##################################################################################################################
+### Plain vanilla estimation
+
+# functional of the distribution to be estimated
+G_fX = a *( m_Y^2 + s_Y^2 - m_Y) + ( 1 - a ) * ( exp( 2 * m_Z + 2 * s_Z ^2) - exp( m_Z + 0.5 * s_Z^2 ));
+print(G_fX);
+
+# series generated by "nature": do not know the distribution
+P = runif(T);
+i_T = t( matrix (QuantileMixture( P, a, m_Y, s_Y, m_Z, s_Z ) ) );
+
+G_Hat_a = function(X) (X[ , 1] - X[ , ncol(X) ]) * X[ , 2 ] * X[ , 2 ];
+G_Hat_b = function(X) apply( X, 1, mean);
+G_Hat_c = function(X) 5 + 0 * X[ , 1];
+G_Hat_d = function(X) (dim(X)[1]-1)/dim(X)[1] * apply( X, 1, var) + apply(X, 1, mean)^2 - apply( X, 1, mean);
+
+Ga = G_Hat_a(i_T); # tentative estimator of unknown functional
+Gb = G_Hat_b(i_T); # tentative estimator of unknown functional
+Gc = G_Hat_c(i_T); # tentative estimator of unknown functional
+Gd = G_Hat_d(i_T); # tentative estimator of unknown functional
+print(Ga);
+print(Gb);
+print(Gc);
+print(Gd);
+
+##################################################################################################################
+### replicability vs. "luck"
+
+# functional of the distribution to be estimated
+G_fX = a *( m_Y^2 + s_Y^2 - m_Y ) + ( 1 - a ) * ( exp( 2 * m_Z + 2*s_Z^2 ) - exp( m_Z + 0.5 * s_Z^2 ) );
+
+# randomize series generated by "nature" to check replicability
+nSim = 10000;
+I_T = matrix( NaN, nSim, T);
+for( t in 1 : T )
+{
+ P = matrix( runif(nSim), nSim, 1);
+ I_T[ , t ] = QuantileMixture( P, a, m_Y, s_Y, m_Z, s_Z );
+}
+
+Ga = G_Hat_a(I_T); # tentative estimator of unknown functional
+Gb = G_Hat_b(I_T); # tentative estimator of unknown functional
+Gc = G_Hat_c(I_T); # tentative estimator of unknown functional
+Gd = G_Hat_d(I_T); # tentative estimator of unknown functional
+
+Loss_Ga = (Ga-G_fX)^2;
+Loss_Gb = (Gb-G_fX)^2;
+Loss_Gc = (Gc-G_fX)^2;
+Loss_Gd = (Gd-G_fX)^2;
+
+Err_Ga = sqrt(mean(Loss_Ga));
+Err_Gb = sqrt(mean(Loss_Gb));
+Err_Gc = sqrt(mean(Loss_Gc));
+Err_Gd = sqrt(mean(Loss_Gd));
+
+Bias_Ga = abs(mean(Ga)-G_fX);
+Bias_Gb = abs(mean(Gb)-G_fX);
+Bias_Gc = abs(mean(Gc)-G_fX);
+Bias_Gd = abs(mean(Gd)-G_fX);
+
+Ineff_Ga = sd(Ga);
+Ineff_Gb = sd(Gb);
+Ineff_Gc = sd(Gc);
+Ineff_Gd = sd(Gd);
+
+##################################################################################################################
+dev.new();
+NumBins = round(10 * log(nSim));
+
+par( mfrow = c(4, 1) );
+
+hist( Ga, NumBins, main = "estimator a" );
+points( G_fX, 0, pch = 21, bg = "red" );
+
+hist( Gb, NumBins, main = "estimator b", xlim = c(-0.2, 1.2) );
+points( G_fX, 0, pch = 21, bg = "red" );
+
+hist( Gc, NumBins, main = "estimator c", xlim = c(-60,60) );
+points( G_fX, 0, pch = 21, bg = "red" );
+
+hist( Gd, NumBins, main = "estimator d" );
+points( G_fX, 0, pch = 21, bg = "red" );
+
+
+#loss
+dev.new();
+par( mfrow = c(4, 1) );
+
+hist( Loss_Ga, NumBins, main = "Loss of estimator a" );
+
+hist( Loss_Gb, NumBins, main = "Loss of estimator b" );
+
+hist( Loss_Gc, NumBins, main = "Loss of estimator c", xlim = c(-60,60) );
+
+hist( Loss_Gd, NumBins, main = "Loss of estimator d" );
+
+##################################################################################################################
+### Stress test replicability
+m_s = seq( 0, 0.2, 0.02 );
+
+Err_Gasq = NULL; Bias_Gasq = NULL; Ineff_Gasq = NULL;
+Err_Gbsq = NULL; Bias_Gbsq = NULL; Ineff_Gbsq = NULL;
+Err_Gcsq = NULL; Bias_Gcsq = NULL; Ineff_Gcsq = NULL;
+Err_Gdsq = NULL; Bias_Gdsq = NULL; Ineff_Gdsq = NULL;
+
+for( j in 1 : length(m_s) )
+{
+ m_Y = m_s[ j ];
+ # functional of the distribution to be estimated
+ G_fX = a * ( m_Y ^ 2 + s_Y^2 - m_Y ) + ( 1 - a ) *( exp( 2 * m_Z + 2 * s_Z^2 ) - exp( m_Z + 0.5 * s_Z^2 ) );
+ # randomize series generated by "nature" to check replicability
+ nSim = 10000;
+ I_T = matrix( NaN, nSim, T);
+ for( t in 1 : T )
+ {
+ P = matrix( runif(nSim) );
+ I_T[ , t ] = QuantileMixture(P,a,m_Y,s_Y,m_Z,s_Z);
+ }
+
+ Ga = G_Hat_a(I_T);
+ Gb = G_Hat_b(I_T);
+ Gc = G_Hat_c(I_T);
+ Gd = G_Hat_d(I_T);
+
+ Loss_Ga = (Ga-G_fX)^2;
+ Loss_Gb = (Gb-G_fX)^2;
+ Loss_Gc = (Gc-G_fX)^2;
+ Loss_Gd = (Gd-G_fX)^2;
+
+ Err_Ga = sqrt(mean(Loss_Ga));
+ Err_Gb = sqrt(mean(Loss_Gb));
+ Err_Gc = sqrt(mean(Loss_Gc));
+ Err_Gd = sqrt(mean(Loss_Gd));
+
+ Bias_Ga = abs(mean(Ga)-G_fX);
+ Bias_Gb = abs(mean(Gb)-G_fX);
+ Bias_Gc = abs(mean(Gc)-G_fX);
+ Bias_Gd = abs(mean(Gd)-G_fX);
+
+ Ineff_Ga = sd(Ga);
+ Ineff_Gb = sd(Gb);
+ Ineff_Gc = sd(Gc);
+ Ineff_Gd = sd(Gd);
+
+ #store results
+ Err_Gasq = cbind( Err_Gasq, Err_Ga^2 ); ##ok<*AGROW>
+ Err_Gbsq = cbind( Err_Gbsq, Err_Gb^2 );
+ Err_Gcsq = cbind( Err_Gcsq, Err_Gc^2 );
+ Err_Gdsq = cbind( Err_Gdsq, Err_Gd^2 );
+
+ Bias_Gasq = cbind(Bias_Gasq, Bias_Ga^2 );
+ Bias_Gbsq = cbind(Bias_Gbsq, Bias_Gb^2 );
+ Bias_Gcsq = cbind(Bias_Gcsq, Bias_Gc^2 );
+ Bias_Gdsq = cbind(Bias_Gdsq, Bias_Gd^2 );
+
+ Ineff_Gasq = cbind( Ineff_Gasq, Ineff_Ga^2 );
+ Ineff_Gbsq = cbind( Ineff_Gbsq, Ineff_Gb^2 );
+ Ineff_Gcsq = cbind( Ineff_Gcsq, Ineff_Gc^2 );
+ Ineff_Gdsq = cbind( Ineff_Gdsq, Ineff_Gd^2 );
+}
+
+##################################################################################################################
+dev.new();
+par( mfrow = c(4, 1) );
+
+
+
+b = barplot( Bias_Gasq + Ineff_Gasq, col = "red", main = "stress-test of estimator a" );
+barplot( Ineff_Gasq, col = "blue", add = TRUE);
+lines( b, Err_Gasq);
+legend( "topleft", 1.9, c( "bias²", "ineff²", "error²" ), col = c( "red","blue", "black" ),
+ lty=1, lwd=c(5,5,1),bg = "gray90" );
+
+b = barplot( Bias_Gbsq + Ineff_Gbsq, col = "red", main = "stress-test of estimator b" );
+barplot( Ineff_Gbsq, col = "blue", add = TRUE);
+lines( b, Err_Gbsq);
+legend( "topleft", 1.9, c( "bias²", "ineff²", "error²" ), col = c( "red","blue", "black" ),
+ lty=1, lwd=c(5,5,1),bg = "gray90" );
+
+b = barplot( Bias_Gcsq + Ineff_Gcsq, col = "red", main = "stress-test of estimator c" );
+barplot( Ineff_Gcsq, col = "blue", add = TRUE);
+lines( b, Err_Gcsq);
+legend( "topleft", 1.9, c( "bias²", "ineff²", "error²" ), col = c( "red","blue", "black" ),
+ lty=1, lwd=c(5,5,1),bg = "gray90" );
+
+b = barplot( Bias_Gdsq + Ineff_Gdsq, col = "red", main = "stress-test of estimator d" );
+barplot( Ineff_Gdsq, col = "blue", add = TRUE);
+lines( b, Err_Gdsq);
+legend( "topleft", 1.9, c( "bias²", "ineff²", "error²" ), col = c( "red","blue", "black" ),
+ lty=1, lwd=c(5,5,1),bg = "gray90" );
\ No newline at end of file
Added: pkg/Meucci/demo/S_Toeplitz.R
===================================================================
--- pkg/Meucci/demo/S_Toeplitz.R (rev 0)
+++ pkg/Meucci/demo/S_Toeplitz.R 2013-07-29 12:18:31 UTC (rev 2662)
@@ -0,0 +1,35 @@
+#' This script shows that the eigenvectors of a Toeplitz matrix have a Fourier basis structure under t-distribution
+#' assumptions, 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_Toeplitz.R"
+#'
+#' @author Xavier Valls \email{flamejat@@gmail.com}
+
+
+###############################################################################################################
+### Inputs
+
+N = 200; # dimension of the matrix
+Decay = 0.9; # decay factor
+
+###############################################################################################################
+T = diag( 1, N);
+for( n in 1 : (N - 1) )
+{
+
+ T = T + Decay^n * ( cbind( matrix( 0, N, N -( N-n ) ), diag( 1, N , N-n) ) +
+ cbind( rbind( matrix(0, N-(N-n), N-n ), diag( 1, N-n)), matrix(0, N, N-(N-n) ) )) ;
+
+
+}
+eig = eigen( T );
+
+###############################################################################################################
+
+#R sorts the eigen vectors, so the results aren't going to be exactly the same as in MATLAB
+
+dev.new();
+plot( eig$vectors[ , n ], type = "l", col = runif(1)*100 );
+lines( eig$vectors[ , n-1 ], type = "l", col = runif(1)*100 );
Added: pkg/Meucci/demo/S_VolatilityClustering.R
===================================================================
--- pkg/Meucci/demo/S_VolatilityClustering.R (rev 0)
+++ pkg/Meucci/demo/S_VolatilityClustering.R 2013-07-29 12:18:31 UTC (rev 2662)
@@ -0,0 +1,30 @@
+#' This file generates paths for a volatility clustering, 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_VolatilityClustering.m"
+#'
+#' @author Xavier Valls \email{flamejat@@gmail.com}
+##################################################################################################################
+### Input parameters
+mu = 0.05; # mean
+a = 0.03;
+b = 0.96;
+s = 0.01;
+T = 1000;
+
+##################################################################################################################
+### Simulate path
+z = rnorm(T);
+s2 = s^2;
+eps = array( NaN, T );
+eps[ 1 ] = s2;
+for( t in 1 : (T - 1) )
+{
+ s2[ t + 1 ] = s^2 + a * ( z[ t ]^2) + b * s2[ t ];
+ eps[ t + 1 ] = mu + sqrt( s2[ t + 1] ) * z[ t + 1 ];
+}
+
+dev.new();
+plot(eps, type = "l", main = "GARCH(1,1) process vs. time", xlab = "", ylab = "" );
Modified: pkg/Meucci/man/Central2Raw.Rd
===================================================================
--- pkg/Meucci/man/Central2Raw.Rd 2013-07-29 10:37:59 UTC (rev 2661)
+++ pkg/Meucci/man/Central2Raw.Rd 2013-07-29 12:18:31 UTC (rev 2662)
@@ -3,18 +3,12 @@
\title{Transforms first n central moments into first n raw moments (first central moment defined as expectation)}
\usage{
Central2Raw(mu)
-
- Central2Raw(mu)
}
\arguments{
- \item{mu}{a vector of central moments}
-
\item{mu}{: [vector] (length N corresponding to order N)
central moments}
}
\value{
- mu_ a vector of non-central moments
-
mu_ : [vector] (length N corresponding to order N)
corresponding raw moments
}
@@ -23,8 +17,6 @@
step 1, we compute the non-central moments. To do so we
start with the first non-central moment and apply
recursively an identity (formula 20)
-
- Map central moments into raw moments
}
\details{
\deqn{ \tilde{ \mu }^{ \big(1\big) }_{X} \equiv \mu
@@ -34,14 +26,13 @@
}
\author{
Ram Ahluwalia \email{rahluwalia at gmail.com}
-
- Xavier Valls \email{flamejat at gmail.com}
}
\references{
A. Meucci - "Exercises in Advanced Risk and Portfolio
Management". See page 10. Symmys site containing original
MATLAB source code \url{http://www.symmys.com}
- \url{http://} See Meucci's script for "Central2Raw.m"
+ \url{http://symmys.com/node/170} See Meucci's script for
+ "Central2Raw.m"
}
Modified: pkg/Meucci/man/Cumul2Raw.Rd
===================================================================
--- pkg/Meucci/man/Cumul2Raw.Rd 2013-07-29 10:37:59 UTC (rev 2661)
+++ pkg/Meucci/man/Cumul2Raw.Rd 2013-07-29 12:18:31 UTC (rev 2662)
@@ -3,28 +3,17 @@
\title{Map cumulative moments into raw moments.}
\usage{
Cumul2Raw(ka)
-
- Cumul2Raw(ka)
}
\arguments{
\item{ka}{: [vector] (length N corresponding to order N)
cumulative moments}
-
- \item{ka}{: [vector] (length N corresponding to order N)
- cumulative moments}
}
\value{
mu_ : [vector] (length N corresponding to order N)
corresponding raw moments
-
- mu_ : [vector] (length N corresponding to order N)
- corresponding raw moments
}
\description{
step 5 of the projection process:
-
- Map cumulative moments into raw moments, as described in
- A. Meucci "Risk and Asset Allocation", Springer, 2005
}
\details{
From the cumulants of Y we compute the raw non-central
@@ -38,10 +27,7 @@
\kappa_{Y}^{ \big(k\big) } \tilde{ \mu } ^{n-k}_{Y} }
}
\author{
- Xavier Valls \email{flamejat at gmail.com} and Ram Ahluwalia
- \email{rahluwalia at gmail.com}
-
- Xavier Valls \email{flamejat at gmail.com}
+ Ram Ahluwalia \email{rahluwalia at gmail.com}
}
\references{
\url{http://symmys.com/node/170} See Meucci's script for
@@ -51,8 +37,5 @@
Skewness, Kurtosis and All Summary Statistics" - formula
(24) Symmys site containing original MATLAB source code
\url{http://www.symmys.com/node/136}
-
- \url{http://symmys.com/node/170} See Meucci's script for
- "Cumul2Raw.m"
}
Modified: pkg/Meucci/man/Raw2Central.Rd
===================================================================
--- pkg/Meucci/man/Raw2Central.Rd 2013-07-29 10:37:59 UTC (rev 2661)
+++ pkg/Meucci/man/Raw2Central.Rd 2013-07-29 12:18:31 UTC (rev 2662)
@@ -3,27 +3,17 @@
\title{Transforms the first n raw moments into the first n central moments}
\usage{
Raw2Central(mu_)
-
- Raw2Central(mu_)
}
\arguments{
- \item{mu_}{the raw (multi-period) non-central moment of
- Y-t}
-
\item{mu_}{: [vector] (length N corresponding to order N)
corresponding raw moments}
}
\value{
- mu (multi-period) central moment of Y-t
-
mu : [vector] (length N corresponding to order N) central
moments
}
\description{
step 6 of projection process:
-
- Map raw moments into central moments, as described in A.
- Meucci "Risk and Asset Allocation", Springer, 2005
}
\details{
compute multi-period central moments.
@@ -37,8 +27,6 @@
}
\author{
Ram Ahluwalia \email{rahluwalia at gmail.com}
-
- Xavier Valls \email{flamejat at gmail.com}
}
\references{
A. Meucci - "Exercises in Advanced Risk and Portfolio
Modified: pkg/Meucci/man/Raw2Cumul.Rd
===================================================================
--- pkg/Meucci/man/Raw2Cumul.Rd 2013-07-29 10:37:59 UTC (rev 2661)
+++ pkg/Meucci/man/Raw2Cumul.Rd 2013-07-29 12:18:31 UTC (rev 2662)
@@ -3,27 +3,18 @@
\title{Transforms raw moments into cumulants}
\usage{
Raw2Cumul(mu_)
-
- Raw2Cumul(mu_)
}
\arguments{
- \item{mu_}{non-central moments of the invariant X-t}
-
\item{mu_}{: [vector] (length N corresponding to order N)
corresponding raw moments}
}
\value{
- ka cumulants of X-t
-
ka : [vector] (length N corresponding to order N)
cumulative moments
}
\description{
Step 3 of the projection process: From the non-central
moments of X-t, we compute the cumulants.
-
- Map raw moments into cumulative moments, as described in
- A. Meucci "Risk and Asset Allocation", Springer, 2005
}
\details{
This process follows from the Taylor approximations for
@@ -38,8 +29,6 @@
}
\author{
Ram Ahluwalia \email{rahluwalia at gmail.com}
-
- Xavier Valls \email{flamejat at gmail.com}
}
\references{
A. Meucci - "Annualization and General Projection of
More information about the Returnanalytics-commits
mailing list