demo/S_EstimateQuantileEvaluation.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 183 - Moment-based functional of a mixture IV".
#'
#' See Meucci's script for "S_EstimateQuantileEvaluation.m"
#'
#' @author Xavier Valls \email{flamejat@@gmail.com}

##################################################################################################################
### Inputs

T   = 52;
a   = 0.5;
m_Y = 0.1;
s_Y = 0.2;
m_Z = 0;
s_Z = 0.15;

##################################################################################################################
### plain estimation

# functional of the distribution to be estimated
G_fX = QuantileMixture( 0.5, a, m_Y, s_Y, m_Z, s_Z);
print( G_fX );

# series generated by "nature": do not know the distribution
P   = runif( T );
i_T = QuantileMixture( P, a, m_Y, s_Y, m_Z, s_Z );

G_Hat_e = function(X) apply( X, 1, median );
G_Hat_b = function(X) apply( X, 1, mean );

Ge = G_Hat_e( t( as.matrix( i_T )) );        # tentative estimator of unknown functional
Gb = G_Hat_b( t( as.matrix( i_T )) );        # tentative estimator of unknown functional
print(Ge);
print(Gb);

##################################################################################################################
### replicability vs. "luck"

# functional of the distribution to be estimated
G_fX = QuantileMixture( 0.5, a, m_Y, s_Y, m_Z, s_Z);

# randomize series generated by "nature" to check replicability
nSim = 10000;
I_T  = c();

for( t in 1:T )
{
    P     = runif(nSim);
    Simul = QuantileMixture(P,a,m_Y,s_Y,m_Z,s_Z);
    I_T   = cbind( I_T, Simul );
}

Ge = G_Hat_e( I_T );        # tentative estimator of unknown functional
Gb = G_Hat_b( I_T );        # tentative estimator of unknown functional

Loss_Ge  = ( Ge - G_fX ) ^ 2;
Loss_Gb  = ( Gb - G_fX ) ^ 2;

Err_Ge   = sqrt( mean( Loss_Ge));
Err_Gb   = sqrt( mean( Loss_Gb));

Bias_Ge  = abs(mean(Ge)-G_fX);
Bias_Gb  = abs(mean(Gb)-G_fX);

Ineff_Ge = sd(Ge);
Ineff_Gb = sd(Gb);

###################################################################################################################
dev.new();
NumBins = round( 10 * log( nSim ));
par( mfrow = c( 2, 1) );
hist( Ge, NumBins, main = "estimator e" );
points(G_fX, 0, pch = 21, bg = "red" );

hist(Gb, NumBins, main = "estimator b" );
points(G_fX, 0, pch = 21, bg = "red" );

#loss
dev.new();
par( mfrow = c( 2, 1) );
hist(Loss_Ge, NumBins, main = "loss of estimator e" );
hist(Loss_Gb, NumBins, main = "loss of estimator b" );

###################################################################################################################
### stress test replicability
m_s = seq( 0, 0.2, 0.02 );

Err_Gesq = NULL; Bias_Gesq = NULL; Ineff_Gesq = NULL;
Err_Gbsq = NULL; Bias_Gbsq = NULL; Ineff_Gbsq = NULL;

for( j in 1 : length( m_s ) )
{
    m_Y = m_s[ j ];
    # functional of the distribution to be estimated
    G_fX = QuantileMixture( 0.5, a, m_Y, s_Y, m_Z, s_Z );

    # randomize series generated by "nature" to check replicability
    nSim = 10000;
    I_T  = NULL;
    for( t in 1 : T )
    {
        P     = runif(nSim);
        Simul = QuantileMixture(P, a, m_Y, s_Y, m_Z, s_Z);
        I_T   = cbind( I_T, Simul );
    }

    Ge = G_Hat_e( I_T );        # tentative estimator of unknown functional
    Gb = G_Hat_b( I_T );        # tentative estimator of unknown functional

    Loss_Ge  = ( Ge - G_fX ) ^ 2;
    Loss_Gb  = ( Gb - G_fX ) ^ 2;

    Err_Ge   = sqrt( mean( Loss_Ge ) );
    Err_Gb   = sqrt( mean( Loss_Gb ) );

    Bias_Ge  = abs( mean( Ge ) - G_fX );
    Bias_Gb  = abs( mean( Gb ) - G_fX );

    Ineff_Ge = std( Ge );
    Ineff_Gb = std( Gb );

    #store results
    Err_Gesq   = cbind( Err_Gesq, Err_Ge^2 ); 
    Err_Gbsq   = cbind( Err_Gbsq, Err_Gb^2 );
    
    Bias_Gesq  = cbind( Bias_Gesq, Bias_Ge^2 ); 
    Bias_Gbsq  = cbind( Bias_Gbsq, Bias_Gb^2 ); 
    
    Ineff_Gesq = cbind( Ineff_Gesq, Ineff_Ge^2 ); 
    Ineff_Gbsq = cbind( Ineff_Gbsq, Ineff_Gb^2 ); 
}

###################################################################################################################
### printlay results
dev.new();
par( mfrow = c( 2, 1) );
b = barplot(Bias_Gesq +Ineff_Gesq , col = "red", main = "stress-test of estimator e");
barplot( Ineff_Gesq, col="blue", add = TRUE);
lines( b, Err_Gesq);
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" );