% =======================================================================
% This is the main run file for estimating the UCSV model for inflation in
% Knüppel and Krüger (JAE, 2021). The codes are based on the codes of
% Joshua Chan used in his paper
% Chan, J.C.C. (2013). Moving Average Stochastic Volatility Models 
% with Application to Inflation Forecast, Journal of Econometrics, 
% 176 (2), 162-172.
% Some more details on the modifications can be found in the file UCSV_noMA.m
% =======================================================================


clear; clc;
if not(isfolder('ucsv_forecasts'))
    mkdir('ucsv_forecasts');
end
addpath('.\matlab_functions');

export_flag = 1;
export_draws_flag = 1;
export_additional_results_flag = 0;

nloop = 15000;
burnin = 5000;

t_start = 60 ;% used for paper
horizmax = 6; % used for paper, but only results for 5 horizons reported there

tic
lags_in_model = 2;
%% load data
% inflation = xlsread('GDPDEF',1,'d16:d297');
price_levels = xlsread('PQvQd.xlsx',1,'b2:hk293');
[~,cols_price_levels] = size(price_levels);
inflation_all = 400*log(price_levels(2:end,:)./price_levels(1:end-1,:));
[~,indx_1st_nan_in_1st_publication] = max(isnan(inflation_all(:,1)));
for i = 2:cols_price_levels 
    inflation_1st_vintage(indx_1st_nan_in_1st_publication+i-2,1) = inflation_all(indx_1st_nan_in_1st_publication+i-2,i); 
    if isnan(inflation_1st_vintage(indx_1st_nan_in_1st_publication+i-2,1))
        inflation_1st_vintage(indx_1st_nan_in_1st_publication+i-2,1) = inflation_all(indx_1st_nan_in_1st_publication+i-2,i+1); %take "second" vintage, only necessary in one case - non-publication in 96Q1 for 95Q4
    end 
end
inflation_1st_vintage(~isnan(inflation_all(:,1))) = inflation_all(~isnan(inflation_all(:,1)),1);
inflation = inflation_1st_vintage(4:285,1); % Use values from 48Q1 to 2018Q2

pfad = '.\ucsv_forecasts\';

t_end = size(inflation,1) ;
meanf = NaN(t_end-2,horizmax);
varf =  NaN(t_end-2,horizmax);
meanf_u = NaN(t_end-2,horizmax);
varf_u =  NaN(t_end-2,horizmax);

y_fore_all = NaN( nloop-burnin , (t_end-t_start+1), horizmax );
mean_y_fore_all = NaN( nloop-burnin , (t_end-t_start+1), horizmax );
vari_y_fore_all = NaN( nloop-burnin , (t_end-t_start+1), horizmax );

tid_all = ((1948 + lags_in_model*0.25):0.25:2018.25)';
counter_tt = 0;

for tt = t_start:t_end   
    counter_tt = counter_tt + 1;
    
    y = inflation(lags_in_model+1:tt,1); 
    T = length(y);
    tid = tid_all(1:T);
    options = optimset('Display', 'off') ;
    warning off all;
    disp(['Data until ' num2str(tid(T))]);

    y0 = inflation(lags_in_model,1); 
    UCSV_noMA; 
      
    meanf(tt-lags_in_model,:) = mean(y_fore);
    varf(tt-lags_in_model,:) = var(y_fore);
    meanf_u(tt-lags_in_model,:) = mean(u_fore);
    varf_u(tt-lags_in_model,:) = var(u_fore);

    y_fore_all( : , counter_tt , : ) = y_fore;
    mean_y_fore_all( : , counter_tt , : ) = mean_y_fore;
    vari_y_fore_all( : , counter_tt , : ) = vari_y_fore;
end

if export_draws_flag == 1   
    if export_additional_results_flag == 1
        horizmax_export = horizmax;
    else
        horizmax_export = horizmax-1;
    end
        
    for ihoriz = 1:horizmax_export
        
        mean_y_fore_all_h = squeeze(mean_y_fore_all(:,:,ihoriz));
        vari_y_fore_all_h = squeeze(vari_y_fore_all(:,:,ihoriz));        
        csvwrite([pfad 'mean_y_fore_all_h' num2str(ihoriz) '.csv' ],mean_y_fore_all_h);
        csvwrite([pfad 'vari_y_fore_all_h' num2str(ihoriz) '.csv' ],vari_y_fore_all_h);

        if export_additional_results_flag == 1
            y_fore_all_h = squeeze(y_fore_all(:,:,ihoriz));
            csvwrite([pfad 'y_fore_all_h' num2str(ihoriz) '.csv' ],y_fore_all_h);
            time_vec = tid(t_start-lags_in_model):0.25:tid(end);
            mean1 = mean(y_fore_all_h);
            mean2 = mean(mean_y_fore_all_h);
            var1 = var(y_fore_all_h);
            var2 = mean( (mean_y_fore_all_h - mean(mean_y_fore_all_h)).^2 ) + mean(vari_y_fore_all_h);
            export_mat = [time_vec;mean1;mean2;var1;var2];
            export_names = {'forecast from:';'mean(y_fore)';'mean(means_y_fore)';'var(y_fore)';'var=f(means_y_fore,vari_y_fore)' };

            xlswrite([pfad 'checks' num2str(ihoriz) '.xlsx' ],export_names,1,'a1');
            xlswrite([pfad 'checks' num2str(ihoriz) '.xlsx' ],export_mat  ,1,'b1');
        end
    end
end

if export_additional_results_flag == 1
    fe = forecast_errors(y,meanf);
    meanf_for = forecasts_from2for(meanf);
    varf_for = forecasts_from2for(varf);
    meanf_u_for = forecasts_from2for(meanf_u);
    varf_u_for = forecasts_from2for(varf_u);

    tid_all_export = tid_all(1:(t_end-lags_in_model));   
    xlswrite([pfad 'overview_all_h.xlsx'],{'dates'},1,'a1');
    xlswrite([pfad 'overview_all_h.xlsx'],{'inflation'},1,'b1');
    xlswrite([pfad 'overview_all_h.xlsx'],{'mean forecasts'},1,'c1');
    xlswrite([pfad 'overview_all_h.xlsx'],{'variance forecasts'},1,'i1');
    xlswrite([pfad 'overview_all_h.xlsx'],{'mean forecast errors'},1,'o1');
    xlswrite([pfad 'overview_all_h.xlsx'],[tid_all_export y meanf_for varf_for fe],1,'a2');                          
end

toc