%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%        Moments selected variables;                           %%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

disp(' ')
disp('%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%')
disp('file dwmoments.m')

nsim        = 1000;  % no. of simulations
ZT          = 11;    % No of variables %needs to be equal to size(ldata)
I           = 1;

disp('nsim =')
disp(nsim)

STDEV    = zeros(nsim,ZT);
RelSTDEV_prod   = zeros(nsim,ZT);
RelSTDEV_gdp    = zeros(nsim,ZT);
SKEW            = zeros(nsim,ZT);
KURT            = zeros(nsim,ZT);
corhp_Y         = zeros(nsim,ZT);
corhp_PIE       = zeros(nsim,ZT);

meanqoq=zeros(nsim,ZT);
medianqoq=zeros(nsim,ZT);
stdqoq=zeros(nsim,ZT);
skewqoq=zeros(nsim,ZT);
kurtqoq=zeros(nsim,ZT);

meanyoy=zeros(nsim,ZT);
medianyoy=zeros(nsim,ZT);
stdyoy=zeros(nsim,ZT);
skewyoy=zeros(nsim,ZT);
kurtyoy=zeros(nsim,ZT);

for sim = 1:nsim;
    options_.simul_seed = 1000+sim*1100;
    options_.periods = 566; %%%NUMBER OF PERIODS to have dif starting points
    options_.drop       = 400;
    stoch_simul(M_,options_,oo_, var_list_); %%%INJI AMMENDED was initially 
    % only var_list_ %%%

clear dataraw;
clear ldata;
clear logdata;

%Transform simulation output into trending data
growth   = g*ones(options_.periods,1);
lPROD    = cumprod(growth);
lOUT     = lPROD.*ytot;
lCONS    = lPROD.*c;
lINV     = lPROD.*inv;
lLABPROD = lPROD.*ytot./n;
lWR      = lPROD.*wr;

ldata        = [lLABPROD lOUT n ur pie piew lWR nint lINV v v./u];
ldata        = ldata(options_.drop+M_.maximum_lag+1:end,:);
ldatanames    = strvcat('labPro','lOUT','n','ur','pie','piew','lWR','nint','linv','v','tight');
[t V]        = size(ldata);


%% Moments HP detrended 
logdata      = log(ldata);

trenddata    = hpfilter(logdata, 1600);
cycledata    = logdata-trenddata;

%HP detrended
STDEV(sim,:)           = nanstd(cycledata,0)*100; % volatilities in percentage deviations from HP
RelSTDEV_prod(sim,:)   = STDEV(sim,:)./STDEV(sim,1);
RelSTDEV_gdp(sim,:)    = STDEV(sim,:)./STDEV(sim,2);
SKEW(sim,:)            = skewness(cycledata,0);
KURT(sim,:)            = kurtosis(cycledata,0)-3;

% Computing correlations between variables:
corhp=[];
    corhp=corrcoef(cycledata(:,:));
clear time;

% Correlation of all variables with GDP 
    for v=1:V
    corhp_Y(sim,v)=corhp(v,2);
    end

% Correlation of all variables with pie  
    for v=1:V
    corhp_PIE(sim,v)=corhp(v,5);
    end
      

%% Moments of qoq and yoy growth rates

dataqoq = 100*[nan(1,size(ldata,2)); ldata(2:end,:)./ldata(1:end-1,:) - 1];
piepos  = strmatch('pie',ldatanames,'exact');
piewpos = strmatch('piew',ldatanames,'exact');
vpos    = strmatch('v',ldatanames,'exact');
dataqoq(2:end,[vpos]) = 100*ldata(2:end,[vpos]); % v are in levels
dataqoq(2:end,[piepos piewpos]) = 100*(ldata(2:end,[piepos piewpos])-1); % price and wage inflation

datayoy = 100*[nan(4,size(ldata,2)); ldata(5:end,:)./ldata(1:end-4,:) - 1]; 
datatmp = nan(size(ldata));
for i = 5:length(ldata)
    datatmp(i,piepos)  = prod(ldata(i-3:i,piepos ))-1; %price inflation
    datatmp(i,piewpos) = prod(ldata(i-3:i,piewpos))-1; %wage inflation
    datatmp(i,vpos)    = ldata(i,vpos)+ldata(i-1,vpos)+ldata(i-2,vpos)+ldata(i-3,vpos); % v are in levels
end
datayoy(5:end,[piepos piewpos vpos]) = 100*datatmp(5:end,[piepos piewpos vpos]);

meanqoq(sim,:) =nanmean(dataqoq);
medianqoq(sim,:) = nanmedian(dataqoq);
stdqoq(sim,:)=nanstd(dataqoq);
skewqoq(sim,:)=skewness(dataqoq);
kurtqoq(sim,:)=kurtosis(dataqoq)-3;

datayoy=datayoy(5:end,:);

meanyoy(sim,:) =nanmean(datayoy);
medianyoy(sim,:) = nanmedian(datayoy);
stdyoy(sim,:)=nanstd(datayoy);
skewyoy(sim,:)=skewness(datayoy);
kurtyoy(sim,:)=kurtosis(datayoy)-3;

end

%% Tables Moments HP-detrended

momdataraw      = [mean(STDEV); mean(RelSTDEV_prod); mean(RelSTDEV_gdp); mean(SKEW); mean(KURT)];
%%%%%%%%%
disp(' ')
%%%%%%%%%
title1=['Model Results - in % deviations from HP'];
labels =strvcat('StDev','RelStDev_prod','RelStDev_gdp','Skewness','Kurtosis');
table(title1,strvcat('var',ldatanames),labels,momdataraw,size(ldatanames,1)+2,2,4);

title2=['Correlation with GDP'];
corY=mean(corhp_Y);
label2 = strvcat('Corr y');
table(title2, strvcat('var',ldatanames),label2,corY,size(ldatanames,1)+2,2,4)

corPIE=mean(corhp_PIE);
title3=['Correlation with PIE'];
label3 = strvcat('Corr pie');
table(title3, strvcat('var',ldatanames),label3,corPIE,size(ldatanames,1)+2,2,4)

%% Tables Moments in growth rates

momQOQ      = [mean(meanqoq); mean(medianqoq);mean(stdqoq); mean(skewqoq);mean(kurtqoq)];
momYOY      =  [mean(meanyoy); mean(medianyoy);mean(stdyoy); mean(skewyoy); mean(kurtyoy)];

labels =strvcat('Mean','Median','StDev','skewness','kurtosis');

title5=['moments of qoq growth rates'];
table(title5,strvcat('var',ldatanames),labels,momQOQ,size(ldatanames,2)+0,0,4);

title6=['moments of yoy growth rates'];
table(title6,strvcat('var',ldatanames),labels,momYOY,size(ldatanames,2)+0,0,4);
