%%
% On Complex Economic Scenario Generators: Is Less More?
% by Jean-Francois Begin
%
% This code creates and estimates the four models investigated in my 
% manuscript: 
% - M1: Wilkie-like (1986)
% - M2: Extended Wilkie-like (1995)
% - M3: Ahlgrim, D'Arcy, and Gorvett (2005, ADG)
% - M4: Begin (2021)

startup;
clear all; clc;

%%
%&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
%&                   0. INITIALIZATION OF VARIABLES                       &
%&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&

% The estimation results have been stored and can be loaded if false is
% selected; they can be regenerated if true is selected
RUN_ESTIMATION = true;

% The argument seed deals with the random seed used to generate our
% results
seed = 1;

% Our sample starts in Sept. 1982 and finishes in Jan. 2020 
dates = datenum(1982,9+(1:448),1);

%%
%&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
%&                1. ESTIMATION OF THE WILKIE-LIKE MODEL                  &
%&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&

if RUN_ESTIMATION
  [chains,acceptance,loglikelihood,model] = getChains('Wilkie',dates,seed);
  save('results/chains_Wilkie','chains','acceptance','loglikelihood','model');
  clear chains acceptance loglikelihood model
end

%%
%&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
%&           2. ESTIMATION OF THE EXTENDED WILKIE-LIKE MODEL              &
%&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&

if RUN_ESTIMATION
  [chains,acceptance,loglikelihood,model] = getChains('ExtWilkie',dates,seed);
  save('results/chains_ExtWilkie','chains','acceptance','loglikelihood','model');
  clear chains acceptance loglikelihood model
end

%%
%&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
%&                 3. ESTIMATION OF THE ADG-LIKE MODEL                    &
%&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&

if RUN_ESTIMATION
  [chains,acceptance,loglikelihood,model] = getChains('ADG',dates,seed);
  save('results/chains_ADG','chains','acceptance','loglikelihood','model');
  clear chains acceptance loglikelihood model
end

%%
%&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
%&                 4. ESTIMATION OF THE NEW ESG MODEL                     &
%&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&

if RUN_ESTIMATION
  [chains,acceptance,loglikelihood,model] = getChains('Begin',dates,seed);
  save('results/chains_Begin','chains','acceptance','loglikelihood','model');
  clear chains acceptance loglikelihood model
end

%%
%&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
%&                  5. GENERATE PATHS FOR THE FOUR MODELS                 &
%&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
T = 360;   % 30 years
N = 10000; % In the paper, we use 100000; with 10000 the results are noisy

[pathsW,modelW] = getSimulationWilkie(T,N);
[pathsE,modelE] = getSimulationExtWilkie(T,N);
[pathsA,modelA] = getSimulationADG(T,N);
[pathsB,modelB] = getSimulationBegin(T,N);

%%
%&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
%&                         6. GENERATE A FEW PLOTS                        &
%&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&

set(0,'DefaultFigureColor', 'w');

fig1 = figData(dates);

fig2 = figMonetary(dates);

fig3 = figFactors(dates);

fig4 = figFunnelofDoubtAll(dates,pathsW,pathsE,pathsA,pathsB);

[fig5,fig6,fig7,fig8,fig9,fig10] = figTracePlots();

[fig11] = figExample(dates,pathsB,7);
