S_EstimateExpectedValueEvaluation.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=52; % number of observations in time series

Mu=.1;
Sigma=.2;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% plain estimation

% unknown functional of the distribution to be estimated, in this case the expected value
G_fX = exp(Mu+0.5*Sigma^2)

i_T=lognrnd(Mu,Sigma,1,T);  % series generated by "nature": do not know the distribution

G1=G_Hat_1(i_T)        % estimator of unknown functional G_1=x(1)*x(3)
G2=G_Hat_2(i_T)        % estimator of unknown functional G_1=sample mean

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% replicability vs. "luck"

% unknown functional of the distribution to be estimated, in this case the expected value
G_fX = exp(Mu+0.5*Sigma^2);

Num_Simulations=10000;
I_T=lognrnd(Mu,Sigma,Num_Simulations,T);  % randomize series generated by "nature" to check replicability

G1=G_Hat_1(I_T);        % estimator of unknown functional G_1=x(1)*x(3)
G2=G_Hat_2(I_T);        % estimator of unknown functional G_2=sample mean

Loss_G1=(G1-G_fX).^2;
Loss_G2=(G2-G_fX).^2;

Err_G1=sqrt(mean(Loss_G1));
Err_G2=sqrt(mean(Loss_G2));

Bias_G1=abs(mean(G1)-G_fX);
Bias_G2=abs(mean(G2)-G_fX);

Ineff_G1=std(G1);
Ineff_G2=std(G2);

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

subplot(2,1,1)
hist(G1,NumBins)
hold on
h=plot(G_fX,0,'.');
set(h,'markersize',20,'color','r')
title('estimator: x(1)*x(3)')

subplot(2,1,2)
hist(G2,NumBins)
hold on
h=plot(G_fX,0,'.');
set(h,'markersize',20,'color','r')
title('estimator: sample mean')

%loss
figure
subplot(2,1,1)
hist(Loss_G1,NumBins)
title('loss of estimator: x(1)*x(3)')
subplot(2,1,2)
hist(Loss_G2,NumBins)
title('loss of estimator: sample mean')

close all


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% stress test replicability
Mus=[0: .1 : .7];

Err_G1sq=[];Err_G2sq=[];Bias_G1sq=[];Bias_G2sq=[];Ineff_G1sq=[];Ineff_G2sq=[];
for j=1:length(Mus)
    Mu=Mus(j);

    % unknown functional of the distribution to be estimated, in this case the expected value
    G_fX = exp(Mu+0.5*Sigma^2);

    Num_Simulations=10000;
    I_T=lognrnd(Mu,Sigma,Num_Simulations,T);  % randomize series generated by "nature" to check replicability

    G1=G_Hat_1(I_T);        % estimator of unknown functional G_1=x(1)*x(3)
    G2=G_Hat_2(I_T);        % estimator of unknown functional G_2=sample mean

    Loss_G1=(G1-G_fX).^2;
    Loss_G2=(G2-G_fX).^2;

    Err_G1=sqrt(mean(Loss_G1));
    Err_G2=sqrt(mean(Loss_G2));
    Err_G1sq=[Err_G1sq Err_G1^2]; %store results
    Err_G2sq=[Err_G2sq Err_G2^2];

    Bias_G1=abs(mean(G1)-G_fX);
    Bias_G2=abs(mean(G2)-G_fX);
    Bias_G1sq=[Bias_G1sq Bias_G1^2]; %store results
    Bias_G2sq=[Bias_G2sq Bias_G2^2];

    Ineff_G1=std(G1);
    Ineff_G2=std(G2);
    Ineff_G1sq=[Ineff_G1sq Ineff_G1^2]; %store results
    Ineff_G2sq=[Ineff_G2sq Ineff_G2^2];

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

    subplot(2,1,1)
    hist(G1,NumBins)
    hold on
    h=plot(G_fX,0,'.');
    set(h,'markersize',20,'color','r')
    title('estimator: x(1)*x(3)')

    subplot(2,1,2)
    hist(G2,NumBins)
    hold on
    h=plot(G_fX,0,'.');
    set(h,'markersize',20,'color','r')
    title('estimator: sample mean')

end
figure

subplot(2,1,1)
h=bar(Bias_G1sq'+Ineff_G1sq','r');
hold on
h=bar(Ineff_G1sq');
hold on
h=plot(Err_G1sq);
legend('bias^2','ineff.^2','error^2','location','best')
title('stress-test of estimator: x(1)*x(3)')

subplot(2,1,2)
h=bar(Bias_G2sq'+Ineff_G2sq','r');
hold on
h=bar(Ineff_G2sq');
hold on
h=plot(Err_G2sq);
legend('bias^2','ineff.^2','error^2','location','best')
title('stress-test of estimator sample mean')