demo/S_EstimateMomentsComboEvaluation.R

In R-Finance/Meucci: Collection of functionality ported from the MATLAB code of Attilio Meucci.

#' 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
#' A. Meucci - "Exercises in Advanced Risk and Portfolio Management" \url{http://symmys.com/node/170},
#' "E 181 - Moment-based functional of a mixture II ".
#'
#' 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 ); 
    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^2", "ineff^2", "error^2" ), 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^2", "ineff^2", "error^2" ), 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^2", "ineff^2", "error^2" ), 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^2", "ineff^2", "error^2" ), col = c( "red","blue", "black" ),
     lty=1, lwd=c(5,5,1),bg = "gray90" );