S_EstimateMomentsComboEvaluation.m

% this script familiarizes the user with the evaluation of an estimator:
% replicability, loss, error, bias and inefficiency
% see "Risk and Asset Allocation"- Springer (2005), by A. Meucci

clear;  close all;  clc;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
T=10;

a=.5;
m_Y=.1;
s_Y=.2;
m_Z=0;
s_Z=.15;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% plain 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+.5*s_Z^2)  );

% series generated by "nature": do not know the distribution
P=rand(T,1)';
i_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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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+.5*s_Z^2)  );

% randomize series generated by "nature" to check replicability
Num_Simulations=10000;
I_T=[];
for t=1:T
    P=rand(Num_Simulations,1);
    Simul=QuantileMixture(P,a,m_Y,s_Y,m_Z,s_Z);
    I_T=[I_T Simul];
end

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=std(Ga);
Ineff_Gb=std(Gb);
Ineff_Gc=std(Gc);
Ineff_Gd=std(Gd);

%estimators
figure 
NumBins=round(10*log(Num_Simulations));

subplot(4,1,1)
hist(Ga,NumBins)
hold on
h=plot(G_fX,0,'.');
set(h,'color','r')
title('estimator a')

subplot(4,1,2)
hist(Gb,NumBins)
hold on
h=plot(G_fX,0,'.');
set(h,'color','r')
title('estimator b')

subplot(4,1,3)
hist(Gc,NumBins)
hold on
h=plot(G_fX,0,'.');
set(h,'color','r')
title('estimator c')

subplot(4,1,4)
hist(Gd,NumBins)
hold on
h=plot(G_fX,0,'.');
set(h,'color','r')
title('estimator d')

%loss
figure 
subplot(4,1,1)
hist(Loss_Ga,NumBins)
title('loss of estimator a')
subplot(4,1,2)
hist(Loss_Gb,NumBins)
title('loss of estimator b')
subplot(4,1,3)
hist(Loss_Gc,NumBins)
title('loss of estimator c')
subplot(4,1,4)
hist(Loss_Gd,NumBins)
title('loss of estimator d')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% stress test replicability
m_s=[0 : .02 : 0.2];

Err_Gasq=[]; Bias_Gasq=[]; Ineff_Gasq=[];
Err_Gbsq=[];Bias_Gbsq=[];Ineff_Gbsq=[];
Err_Gcsq=[];Bias_Gcsq=[];Ineff_Gcsq=[];
Err_Gdsq=[];Bias_Gdsq=[];Ineff_Gdsq=[];

for j=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+.5*s_Z^2)  );
    % randomize series generated by "nature" to check replicability
    Num_Simulations=10000;
    I_T=[];
    for t=1:T
        P=rand(Num_Simulations,1);
        Simul=QuantileMixture(P,a,m_Y,s_Y,m_Z,s_Z);
        I_T=[I_T Simul];
    end

    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=std(Ga);
    Ineff_Gb=std(Gb);
    Ineff_Gc=std(Gc);
    Ineff_Gd=std(Gd);

    %store results
    Err_Gasq=[Err_Gasq Err_Ga^2];
    Err_Gbsq=[Err_Gbsq Err_Gb^2];
    Err_Gcsq=[Err_Gcsq Err_Gc^2];
    Err_Gdsq=[Err_Gdsq Err_Gd^2];
    
    Bias_Gasq=[Bias_Gasq Bias_Ga^2]; 
    Bias_Gbsq=[Bias_Gbsq Bias_Gb^2]; 
    Bias_Gcsq=[Bias_Gcsq Bias_Gc^2]; 
    Bias_Gdsq=[Bias_Gdsq Bias_Gd^2]; 
    
    Ineff_Gasq=[Ineff_Gasq Ineff_Ga^2]; 
    Ineff_Gbsq=[Ineff_Gbsq Ineff_Gb^2]; 
    Ineff_Gcsq=[Ineff_Gcsq Ineff_Gc^2]; 
    Ineff_Gdsq=[Ineff_Gdsq Ineff_Gd^2]; 


end

figure

subplot(4,1,1)
h=bar(Bias_Gasq'+Ineff_Gasq','r');
hold on
h=bar(Ineff_Gasq');
hold on
h=plot(Err_Gasq);
legend('bias^2','ineff.^2','error^2','location','best')
title('stress-test of estimator a')

subplot(4,1,2)
h=bar(Bias_Gbsq'+Ineff_Gbsq','r');
hold on
h=bar(Ineff_Gbsq');
hold on
h=plot(Err_Gbsq);
legend('bias^2','ineff.^2','error^2','location','best')
title('stress-test of estimator b')

subplot(4,1,3)
h=bar(Bias_Gcsq'+Ineff_Gcsq','r');
hold on
h=bar(Ineff_Gcsq');
hold on
h=plot(Err_Gcsq);
legend('bias^2','ineff.^2','error^2','location','best')
title('stress-test of estimator c')

subplot(4,1,4)
h=bar(Bias_Gdsq'+Ineff_Gdsq','r');
hold on
h=bar(Ineff_Gdsq');
hold on
h=plot(Err_Gdsq);
legend('bias^2','ineff.^2','error^2','location','best')
title('stress-test of estimator d')