%%%%%%%%%%%
%% FORECAST
%%%%%%%%%%%
close all;
clear all;

addpath('../infrastructure/tools');
addpath('./data');

%% Name baseline
baseline = 'May_MPC';
control  = 'Mar_MPC';

%% Defines the time frame of the forecast
% Don't forget to change the time frame depending on you data and forecast
% period!
startfcast = qq(2018,1);
endMPC     = qq(2020,4);
endfcast   = qq(2022,4);

% forecast range
simrng = startfcast:endfcast;

% judgment range
mpcrng = startfcast:endMPC;

t = length(simrng);
k = length(mpcrng);

%% Load model and filtered historical data
load('../infrastructure/model/model.mat');
load('./data/filter_db.mat');

h = dbclip(filtdb.mean,qq(1990,1):startfcast-1);

%% Load assumptions for hard-tunes over the forecast
a = dbload('data/history.csv');

%% Load core model results for comparison purposes
core     = dbload(['data/', baseline, '_coremodel.csv']); 
prevCore = dbload(['data/', control , '_coremodel_inclVAT.csv']); 
gpm      = dbload(['data/WGDP_Assumption_May 2018 MPC.csv']);
base     = dbload(['data/', baseline, '_base.csv']);
disag    = dbload(['data/', baseline, '_DIM.csv']); 
foodVAR  = dbload(['data/', baseline, '_foodvar.csv']);
prevMPC  = dbload(['data/', control , '_fcastdata.csv']);
highfreq = dbload(['data/', baseline, '_highfreq.csv']);

%% In the forecast simulation variables can be set as numbers (for instance
% for external development as opposed to automatically used AR(1) processes). 
% These preset variables are in the historical database and they require 
% creation of a 'plan' for this variable. The command 'plan' creates plan 
% for model 'm' over the specified horizon in the object 'simplan'. These
% variables must be marked using the command 'exogenize', while the residual
% used for equalizing the model equation for the fixed variable must be marked using
% command 'endogenize'. Note that any variable that is fixed must have residual term
% as a part of its equation. The forecast is simulated in a similar way as 
% without the fixing. See example below where external inflation 
% and external interest rates could be taken, ay, from IMF's WEO or 
% Consensus Forecast. Be aware that you must specify exactly till when
% the data is available in the database.

% In case you use exact numbers for external variables, i.e. foreing inflation,
% interest rate or oil prices, the object "simplan" created below has to be inserted
% into the command simulate.

simplan = plan(m,simrng);

% Specify the exogenize and endogenize commands for vector of external variables
% All (except gdpot6) assumptions are AUTOREGRESSIVE to steady state after jrng

%% ASSUMPTIONS 

% [LCPI_US, LCPI_EU, LCPI_JA]
simplan = exogenize(simplan,{'dot_x_cpi'},simrng);
simplan = endogenize(simplan,{'e_dot_x_cpi'},simrng);
h.dot_x_cpi(mpcrng(1:k)) = a.dot_x_cpi(mpcrng(1:k)); 
h.dot_x_cpi(simrng(k+1:end)) = arf(a.dot_x_cpi,[1 -0.8],...
    m.dot_x_cpi*(1-0.8),simrng(k:end)) ;


% [RS_US, RS_EU, RS_JA]
simplan = exogenize(simplan,{'x_rn'},simrng);
simplan = endogenize(simplan,{'e_x_rn'},simrng);
h.x_rn(mpcrng(1:k)) = a.x_rn(mpcrng(1:k)); 
h.x_rn(simrng(k+1:end)) = arf(a.x_rn,[1 -0.90],...
    m.x_rn*(1-0.90),simrng(k:end)) ;

% Y_WORLD
simplan = exogenize(simplan,{'lx_gdp_gap'},simrng);
simplan = endogenize(simplan,{'e_lx_gdp_gap'},simrng);
h.lx_gdp_gap(mpcrng(1:k)) = a.lx_gdp_gap(mpcrng(1:k));
h.lx_gdp_gap(simrng(k+1:end)) = arf(a.lx_gdp_gap,[1 -0.8],...
    m.lx_gdp_gap*(1-0.8),simrng(k:end)) ;

% WORLD GROWTH
simplan = exogenize(simplan,{'dot_x_gdp'},simrng);
simplan = endogenize(simplan,{'e_lx_gdp_eq'},simrng);
h.dot_x_gdp(mpcrng(1:k)) = a.dot_x_gdp(mpcrng(1:k));
h.dot_x_gdp(simrng(k+1:end)) = arf(a.dot_x_gdp,[1 -0.8],...
    m.dot_x_gdp*(1-0.8),simrng(k+1:end)) ;

% WFOOD08
simplan = exogenize(simplan,{'dot_food'},simrng);
simplan = endogenize(simplan,{'e_dot_qfood_eq'},simrng);
simplan = endogenize(simplan,{'e_lqfood_gap'},simrng);
h.dot_food(mpcrng(1:k)) = a.dot_food(mpcrng(1:k));
h.dot_food(simrng(k+1:end)) = arf(a.dot_food,[1 -0.8],...
    m.dot_food*(1-0.8),simrng(k+1:end)) ;
 
% OILPDD
simplan = exogenize(simplan,{'dot_oil'},simrng);
simplan = endogenize(simplan,{'e_dot_qoil_eq'},simrng);
simplan = endogenize(simplan,{'e_lqoil_gap'},simrng);
h.dot_oil(mpcrng(1:k)) = a.dot_oil(mpcrng(1:k));
h.dot_oil(simrng(k+1:end)) = arf(a.dot_oil,[1 -0.8],...
    m.dot_oil*(1-0.8),simrng(k+1:end)) ;

% PETTAX
simplan = exogenize(simplan,{'dot_pettax'},simrng);
simplan = endogenize(simplan,{'e_dot_pettax'},simrng);
h.dot_pettax(mpcrng(1:k)) = a.dot_pettax(mpcrng(1:k));
h.dot_pettax(simrng(k+1:end)) = arf(a.dot_pettax,[1 -0.9],...
    m.dot_pettax*(1-0.9),simrng(k+1:end)) ;

% CPIELECD08
simplan = exogenize(simplan,{'dot_cpi_elec'},simrng(1:k));
simplan = endogenize(simplan,{'e_dot_cpi_elec'},simrng(1:k));
h.dot_cpi_elec(mpcrng(1:k)) = a.dot_cpi_elec(mpcrng(1:k)); 

simplan = exogenize(simplan,{'dot4_cpi_elec'},simrng(k+1:end));
simplan = endogenize(simplan,{'e_dot_cpi_elec'},simrng(k+1:end));
h.dot4_cpi_elec(simrng(k+1:end)) = 100*log(1+[8]/100);

% GDPOT6 - G
simplan = exogenize(simplan,{'gg'},mpcrng(1:k));
simplan = endogenize(simplan,{'e_g'},mpcrng(1:k));
h.gg(mpcrng(1:k)) = a.dot_gdp_pot(mpcrng(1:k));

% GDPOT6 - G
simplan = exogenize(simplan,{'dot_gdp_eq'},mpcrng(1:k));
simplan = endogenize(simplan,{'e_lgdp_eq'},mpcrng(1:k));
h.dot_gdp_eq(mpcrng(1:k)) = a.dot_gdp_bar(mpcrng(1:k));

% FORCMD
simplan = exogenize(simplan,{'dot_forcmd'},simrng);
simplan = endogenize(simplan,{'e_dot_qforcmd_eq'},simrng);
simplan = endogenize(simplan,{'e_lqforcmd_gap'},simrng);
h.dot_forcmd(mpcrng(1:k)) = a.dot_forcmd(mpcrng(1:k)) ;
h.dot_forcmd(simrng(k+1:end)) = arf(a.dot_forcmd,[1 -0.8],...
    m.dot_forcmd*(1-0.8),simrng(k+1:end)) ;


%% NTF or expert forecast - hard tunes

% Foreign RIR equilibrium fix
% These values are based on a4_driverFilter run to 2020Q4
simplan = exogenize(simplan,{'x_rr_eq'},simrng);
simplan = endogenize(simplan,{'e_x_rr_eq'},simrng);
h.x_rr_eq(simrng) = arf(h.x_rr_eq,[1 -0.95],...
    m.x_rr_eq*(1-0.95),simrng);

% RIR equilibrium fix
simplan = exogenize(simplan,{'rr_eq'},simrng);
simplan = endogenize(simplan,{'e_rr_eq'},simrng);
h.rr_eq(simrng) = arf(h.rr_eq,[1 -0.905],...
    m.rr_eq*(1-0.905),simrng);

% % REER equilibrium fix
simplan = exogenize(simplan,{'dot_z_eq'},simrng(1));
simplan = endogenize(simplan,{'e_dot_z_eq'},simrng(1));
h.dot_z_eq(simrng(1)) = -1.9; 

% FORCMD equilibrium fix
simplan = exogenize(simplan,{'dot_qforcmd_eq'},mpcrng);
simplan = endogenize(simplan,{'e_dot_qforcmd_eq'},mpcrng);
h.dot_qforcmd_eq(mpcrng) = arf(h.dot_qforcmd_eq,[1 -0.6],...
    0.0*(1-0.6),mpcrng);

% [EXDOLLD, LS_EU, LS_JA]
simplan = exogenize(simplan,{'dot_s'}, simrng(1:3));  
simplan = endogenize(simplan,{'e_ls'}, simrng(1:3)); 
h.dot_s(mpcrng(1:3)) = [highfreq.dot_s(qq(2018,1)) highfreq.dot_s(qq(2018,2)) 1.6293]; 

% REPORI
simplan = exogenize(simplan,{'rn'},mpcrng(1:2));
simplan = endogenize(simplan,{'e_rn'},mpcrng(1:2));
h.rn(mpcrng(1:2)) = [6.74 6.5];

% Long-term interest rates
simplan = exogenize(simplan,{'rn8'},mpcrng(1:2));
simplan = endogenize(simplan,{'e_term8'},mpcrng(1:2));
h.rn8(mpcrng(1)) = highfreq.rn8(mpcrng(1));
h.rn8(mpcrng(2)) = highfreq.rn8(mpcrng(2));

simplan = exogenize(simplan,{'rn20'},mpcrng(1:2));
simplan = endogenize(simplan,{'e_term20'},mpcrng(1:2));
h.rn20(mpcrng(1)) = highfreq.rn20(mpcrng(1));
h.rn20(mpcrng(2)) = highfreq.rn20(mpcrng(2));

simplan = exogenize(simplan,{'rn40'},mpcrng(1:2));
simplan = endogenize(simplan,{'e_term40'},mpcrng(1:2));
h.rn40(mpcrng(1)) = highfreq.rn40(mpcrng(1));
h.rn40(mpcrng(2)) = highfreq.rn40(mpcrng(2));

% GDPMP6
simplan = exogenize(simplan,{'dot_gdp'},mpcrng(1:2));
simplan = endogenize(simplan,{'e_lgdp_gap'},mpcrng(1:2));
h.dot_gdp(mpcrng(1:2)) = 100*log(1+[-0.4 2.3]/100);

% SERVICES
simplan = exogenize(simplan,{'dot_cpi_serv'},mpcrng(1:3));
simplan = endogenize(simplan,{'e_dot_cpi_serv'},mpcrng(1:3));
h.dot_cpi_serv(mpcrng(1)) = a.dot_cpi_serv(mpcrng(1));
h.dot_cpi_serv(mpcrng(2)) = prevMPC.dot_cpi_serv(mpcrng(2))-0.4; 
h.dot_cpi_serv(mpcrng(3)) = prevMPC.dot_cpi_serv(mpcrng(3));

% GOODSX
simplan = exogenize(simplan,{'dot_cpi_goodsx'},mpcrng(1:3));
simplan = endogenize(simplan,{'e_dot_cpi_goodsx'},mpcrng(1:3));
h.dot_cpi_goodsx(mpcrng(1)) = a.dot_cpi_goodsx(mpcrng(1)); 
h.dot_cpi_goodsx(mpcrng(2)) = a.dot_cpi_goodsx(mpcrng(2))+0.49; 
h.dot_cpi_goodsx(mpcrng(3)) = a.dot_cpi_goodsx(mpcrng(3))+0.1; 

% CPIFR
simplan = exogenize(simplan,{'dot_cpi_food'},mpcrng(1:2));
simplan = endogenize(simplan,{'e_dot_cpi_food'},mpcrng(1:2));
h.dot_cpi_food(mpcrng(1)) = a.dot_cpi_food(mpcrng(1));
h.dot_cpi_food(mpcrng(2)) = a.dot_cpi_food(mpcrng(2)) + 2.12; %starting point 

% CORE
simplan = exogenize(simplan,{'dot_cpi_core'},mpcrng(1));
simplan = endogenize(simplan,{'e_lcpi_core'},mpcrng(1));
h.dot_cpi_core(mpcrng(1)) = a.dot_cpi_core(mpcrng(1)) ;

% BFP
simplan = exogenize(simplan,{'dot_bfp'},mpcrng(1:3));
simplan = endogenize(simplan,{'e_dot_bfp'},mpcrng(1:3));
h.dot_bfp(mpcrng(1:3)) = a.dot_bfp(mpcrng(1:3));

% CPIPETR 
simplan = exogenize(simplan,{'dot_cpi_petr'},mpcrng(1:3));
simplan = endogenize(simplan,{'e_dot_cpi_petr'},mpcrng(1:3));
h.dot_cpi_petr(mpcrng(1:3)) = a.dot_cpi_petr(mpcrng(1:3));

% EMPLO5
simplan = exogenize(simplan,{'dot_emp'},mpcrng(1));
simplan = endogenize(simplan,{'e_lemp_gap'},mpcrng(1));
h.dot_emp(mpcrng(1)) = 100*log(1+[1.3]/100);

% PREM
simplan = exogenize(simplan,{'prem'},mpcrng(1:2));
simplan = endogenize(simplan,{'e_prem'},mpcrng(1:2));
h.prem(mpcrng(1)) = highfreq.prem(mpcrng(1)); 
h.prem(mpcrng(2)) = highfreq.prem(mpcrng(2)); 

simplan = exogenize(simplan,{'prem_eq'}, simrng);
simplan = endogenize(simplan,{'e_prem_eq'}, simrng);
h.prem_eq(simrng) = arf(h.prem_eq,[1 -0.99],...
     (m.prem_eq)*(1-0.99),simrng);
 
% Weights
simplan = exogenize(simplan,{'w_food'},simrng);
simplan = endogenize(simplan,{'e_w_food'},simrng);
h.w_food(simrng) = 0.1724;

simplan = exogenize(simplan,{'w_petr'},simrng);
simplan = endogenize(simplan,{'e_w_petr'},simrng);
h.w_petr(simrng) = 0.0458;

simplan = exogenize(simplan,{'w_elec'},simrng);
simplan = endogenize(simplan,{'e_w_elec'},simrng);
h.w_elec(simrng) = 0.0375;

simplan = exogenize(simplan,{'w_goodsx'},simrng);
simplan = endogenize(simplan,{'e_w_goodsx'},simrng);
h.w_goodsx(simrng) = 0.2313;

simplan = exogenize(simplan,{'w_serv'},simrng);
simplan = endogenize(simplan,{'e_w_serv'},simrng);
h.w_serv(simrng) = 0.513;

simplan = exogenize(simplan,{'w_bfp'},simrng);
simplan = endogenize(simplan,{'e_w_bfp'},simrng);
h.w_bfp(simrng) = 0.4502;
             
%% Residuals (shocks) - soft tunes

% Automatic - smooth out all residuals that kalman identified on history
% note well that the tuning of residuals for NTF vars starts 
% one period later MAKES THE AR!!!!! enames are the residuals [1 -0.5] are
% the inverse roots [0] is the mean
% nnames = get(m,'enames');
% for ii = 1:length(nnames);
%     h.(nnames{ii}) = arf(h.(nnames{ii}),[1 -0.4],[0],fcastrng);
% end 

% Manual - loaded from shock_plan.csv
shocks=dbload('data\shock_plan.csv');

sh = get(m, 'elist');
for i = 1 : numel(sh)
    h.(sh{i})(mpcrng(1):mpcrng(k)+12) = shocks.(sh{i})(mpcrng(1):mpcrng(k)+12);    
end

%% Switch off Taylor rule modification that ignores direct VAT impacts 
p = get(m,'params');
p.vat_switch = 0;
m = assign(m,p);
m = solve(m);

%% Simulate
% Command 'simulate' simulates the model 'm' based on the database 'h' over
% the forecast range 'fcastrange'. Results are written into the object 's'.
s = simulate(m,h,simrng,'plan=',simplan,'anticipate=',true);
% s1 = jforecast(m,h,simrng,'plan=',simplan);

% Command 'dboverlay' puts together the historical database 'h' and the
% results of the simulation saved in object 's'. Single database 's' is
% created
h = dbextend(h,s);
s = h;

% % Save in-sample shocks to give indication of starting points for the
% % forecast. Paste these historical shocks in shock_plan.csv (only has to be
% % done once, when historical data is final.)
% 
% hist_shocks=[];
% for i = 1 : numel(sh)
%   hist_shocks.(sh{i})= tseries;  
%   hist_shocks.(sh{i})(qq(2015,1):simrng(end))=h.(sh{i})(qq(2015,1):simrng(end));
% end
% dbsave(hist_shocks, 'data/hist_shocks.csv', qq(2016,1):simrng(end));
% 
% return;

%% Calculate actual percentage changes for reporting
s.CPITD       = exp(s.lcpi/100);
s.CPIFD       = exp(s.lcpi_food/100);
s.CPIXFEPD    = exp(s.lcpi_core/100); 
s.CPIPETD     = exp(s.lcpi_petr/100);
s.CPIELECD    = exp(s.lcpi_elec/100);
s.GDPMP6      = exp(s.lgdp/100);
s.GDPOT6      = exp(s.lgdp_eq/100);
s.SALAV5      = exp(s.lwage/100);
s.EMPLO5      = exp(s.lemp/100);
s.ULCSTD      = s.SALAV5*s.EMPLO5/s.GDPMP6;
s.RWAGE       = exp(s.lrw/100);
s.IMPTD       = exp((s.lx_cpi+s.ls)/100);
s.CPISERVD    = exp(s.lcpi_serv/100);
s.CPIGDXD     = exp(s.lcpi_goodsx/100);
s.WGDP        = exp(s.lx_gdp/100);
s.PROD        = exp(s.lgdp/100)/exp(s.lemp/100);
s.BRENT       = exp(s.loil/100);
s.EXREALD     = exp(-s.lz/100);
s.EXNOMD      = exp(-s.ls/100);  

s.CPITR       = ((s.CPITD/s.CPITD{-4})-1)*100;
s.CPIFR       = ((s.CPIFD/s.CPIFD{-4})-1)*100;
s.CPIXFEPR    = ((s.CPIXFEPD/s.CPIXFEPD{-4})-1)*100;
s.CPIPETR     = ((s.CPIPETD/s.CPIPETD{-4})-1)*100;
s.CPIELECR    = ((s.CPIELECD/s.CPIELECD{-4})-1)*100;
s.GDPMPR      = ((s.GDPMP6/s.GDPMP6{-4})-1)*100;
s.GDPOTR      = ((s.GDPOT6/s.GDPOT6{-4})-1)*100;
s.SALAVR      = ((s.SALAV5/s.SALAV5{-4})-1)*100;
s.ULCSTR      = ((s.ULCSTD/s.ULCSTD{-4})-1)*100;
s.EMPLOR      = ((s.EMPLO5/s.EMPLO5{-4})-1)*100;
s.RWAGER      = ((s.RWAGE/s.RWAGE{-4})-1)*100;
s.IMPTR       = ((s.IMPTD/s.IMPTD{-4})-1)*100;
s.CPISERVR    = ((s.CPISERVD/s.CPISERVD{-4})-1)*100;
s.CPIGDXR     = ((s.CPIGDXD/s.CPIGDXD{-4})-1)*100;
s.WGDPR       = ((s.WGDP/s.WGDP{-4})-1)*100;
s.PRODR       = ((s.PROD/s.PROD{-4})-1)*100;
s.BRENTR      = ((s.BRENT/s.BRENT{-4})-1)*100;
s.EXREALDR    = ((s.EXREALD /s.EXREALD{-4})-1)*100;
s.EXNOMDR     = ((s.EXNOMD /s.EXNOMD{-4})-1)*100;

s.CPITRq      = (((s.CPITD/s.CPITD{-1})^4)-1)*100;
s.CPIFRq      = (((s.CPIFD/s.CPIFD{-1})^4)-1)*100;
s.CPIXFEPRq   = (((s.CPIXFEPD/s.CPIXFEPD{-1})^4)-1)*100;
s.GDPMPRq     = (((s.GDPMP6/s.GDPMP6{-1})^4)-1)*100;
s.GDPOTRq     = (((s.GDPOT6/s.GDPOT6{-1})^4)-1)*100;
s.SALAVRq     = (((s.SALAV5/s.SALAV5{-1})^4)-1)*100;
s.RWAGERq     = (((s.RWAGE/s.RWAGE{-1})^4)-1)*100;
s.CPISERVRq   = (((s.CPISERVD/s.CPISERVD{-1})^4)-1)*100;
s.CPIGDXRq    = (((s.CPIGDXD/s.CPIGDXD{-1})^4)-1)*100;
s.WGDPRq      = (((s.WGDP/s.WGDP{-1})^4)-1)*100;
s.EXREALDRq   = (((s.EXREALD /s.EXREALD{-1})^4)-1)*100;
s.EXNOMDRq    = (((s.EXNOMD /s.EXNOMD{-1})^4)-1)*100;


%% Generate Headline CPI for reporting

s.CPITR_SSA = tseries;
s.CPITR_SSA(qq(1990,1):startfcast-1) = core.CPITR(qq(1990,1):startfcast-1);
s.CPITR_SSA(startfcast:endfcast) = s.CPITR(startfcast:endfcast);

s.CPIXFEPR_SSA = tseries;
s.CPIXFEPR_SSA(qq(1990,1):startfcast-1) = core.CPIXFEPR08(qq(1990,1):startfcast-1);
s.CPIXFEPR_SSA(startfcast:endfcast) = s.CPIXFEPR(startfcast:endfcast);

s.CPIPETR_SSA = tseries;
s.CPIPETR_SSA(qq(1990,1):startfcast-1) = core.CPIPETR08(qq(1990,1):startfcast-1);
s.CPIPETR_SSA(startfcast:endfcast) = s.CPIPETR(startfcast:endfcast);

%% Save results

% Results are saved in file 'fcastdata.csv'  
dbsave(['data/', baseline, '_fcastdata.csv'],s);
save(['data/', baseline, '_fcastdata.mat'],'s');

%% Report

% Add zeroline to database
s.zeroline = tseries;
s.zeroline(qq(2000,1):qq(2021,4)) = 0;

% Plot range
prange = qq(2014,1):qq(2020,4);
mpcrng=mpcrng(1:k);

commandwindow;
graphsForecast;

%% Finish

rmpath('../infrastructure/tools');
rmpath('./data');