% adapted from figure3.gss
% compute figure 3


workspacedir='.\workspace\';

%% parameters
nfac=3;                             % Total number of factors %
n_i1_vec=[1;1;0];                     % Indicator for I1 factors %
nar=4;
nma=8;                              % Number of AR lags, including first difference lags %

% Construct Model Identifier String %
id_str=['workspace_' num2str(nfac) '_' num2str(nar) '_' num2str(n_i1_vec)'];

% Number of factors, excluding NPI %
nf=nfac-1;

%% Load Results %
load([workspacedir id_str]);
%loaded: 'em_it_completed','npi_llf','alpha','lam','rho','phi','var_eps','var_e','x3t','p3t','f'
%rename some
llf_sum=npi_llf_sum;
x3t_model=x3t;
p3t_model=p3t;

[nobs,nseries]=size(xdata);
 step=1/4;
 calds=(nfirst:step:nfirst+(nobs-1)*step)';

nstate=size(f,1);
nobs=size(x3t_model,1)-1;
q=zeros(nstate,nstate);
q(1:nfac,1:nfac)=var_e;

%% Transform Model to First differences, with levels added to end of state vector %
g=zeros(nstate,nstate);
for i =1: nar
    j=(i-1)*nfac;
    g(j+1:j+nfac,j+1:j+nfac)=eye(nfac);
    g(j+1:j+nfac,j+nfac+1:j+2*nfac)=-eye(nfac);
end
j=nar*nfac;
g(j+1:j+nfac,j+1:j+nfac)=eye(nfac);
gi=inv(g);
fd=g*f*gi;
qd=g*q*g';
x3t_model=x3t_model*g';
% Delete rows and columns of state that are integrated %
nn=logical([zeros(nfac*nar,1);n_i1_vec]);
f=fd(~nn,~nn);
q=qd(~nn,~nn);
x3t_model=x3t_model(:,~nn');

nur=sum(n_i1_vec);  % Number of Unit Roots %
ntmp=(nstate-nur)*(nstate-nur);
ptmp=zeros(size(p3t_model,1),ntmp);
for t =1: size(x3t_model,1)
    tmp=reshape(p3t_model(t,:),nstate,nstate);
    tmp=g*tmp*g';
    tmp=tmp(~nn,~nn);
    ptmp(t,:)=tmp(:)';
end
p3t_model=ptmp;


%% Compute Variance of a(t) ... from projection of n onto leads and lags
%   of f ... implemented using a Kalman Smoother
%   That is, set up a state space system with the F's observed and
%   compute decomposition of n
% %

nstate=size(f,1);
h=zeros(nstate,nf);
h(2:nf+1,:)=eye(nf);
r=zeros(nf,nf);

% Compute unconditional variance %
np=nstate*nstate;
p=inv(eye(np) - kron(f,f))*q(:);
p=reshape(p,nstate,nstate);

% -- KS -- %
p1t=zeros(nobs+1,np);
p2t=zeros(nobs+1,np);
p3t=zeros(nobs+1,np);
x1t=zeros(nobs+1,nstate);
x2t=zeros(nobs+1,nstate);
x3t=zeros(nobs+1,nstate);
x1=zeros(nstate,1);
p1=p;
x1t(1,:)=x1';
p1t(1,:)=(p1(:))';
for i =1: nobs
    y=x3t_model(i+1,2:nfac)';
    [x1,p1,x2,p2]=kfilt_llf(y,x1,p1,h,f,r,q);
    x1t(i+1,:)=x1';
    p1t(i+1,:)=(p1(:))';
    x2t(i+1,:)=x2';
    p2t(i+1,:)=(p2(:))';
end

% -- Kalman Smoother -- %
x3t(nobs+1,:)=x1';
p3t(nobs+1,:)=(p1(:))';
x3=x1; p3=p1;
for i =nobs:-1:2
    x2=x2t(i+1,:)';
    p2=p2t(i+1,:)';
    p2=reshape(p2,nstate,nstate);
    x1=x1t(i,:)';
    p1=p1t(i,:)';
    p1=reshape(p1,nstate,nstate);
    [x3,p3]=Ksmooth(x1,x2,x3,p1,p2,p3,f);
    x3t(i,:)=x3';
    p3t(i,:)=(p3(:))';
end
x3t_1=x3t;
x1t_1=x1t;

% Construct Decomposition of npi %
dnpi=x3t_model(:,1);
dnpi_fit=x3t(:,1);
dnpi_ex=dnpi-dnpi_fit;

npi=cumsum(dnpi);
npi_fit=cumsum(dnpi_fit);
npi_ex=cumsum(dnpi_ex);

drfac=x3t_model(:,2:nfac);
rfac=cumsum(drfac);
drfac_fit=x3t(:,2:nfac);
rfac_fit=cumsum(drfac_fit);

%tmp=[calds npi npi_fit npi_ex rfac];

npi_ma=nan(size(calds,1),1);
npi_fit_ma=nan(size(calds,1),1);
npi_ex_ma=nan(size(calds,1),1);
rfac_ma=nan(size(calds,1),2);
for i =nma+1: size(dnpi,1)
    npi_ma(i,:)=sum(dnpi(i-nma+1:i,:));
    npi_fit_ma(i,:)=sum(dnpi_fit(i-nma+1:i,:));
    npi_ex_ma(i,:)=sum(dnpi_ex(i-nma+1:i,:));
    rfac_ma(i,:)=sum(drfac(i-nma+1:i,:))';
end


%select PCE data
isel = (cal_mac >= calds(1)-.0001) & (cal_mac < calds(end)+.0001);
y=d_pcep(isel);

x=[4*npi_ex 4*npi_fit rfac];
tmp=[calds y x];
tmp=tmp(~sum(isnan(tmp),2),:);
calp=tmp(:,1);
yp=tmp(:,2);
xp=tmp(:,3:end);

yp=yp-mean(yp);
xp=xp-mean(xp);
rp=xp(:,3:end);
v=xp(:,1);
ar=xp(:,2);
zp=yp-v-ar;
%b=zp/rp;
b=(rp'*rp)\rp'*zp;
ye=zp-rp*b;
yv=v;
yar=ar;
yrp=rp*b;

isel = (calp > yfirst-.000001) & (calp < ylast +.0000001);
calp=calp(isel);
yp=yp(isel);
yv=yv(isel);
yar=yar(isel);
yrp=yrp(isel);
ye=ye(isel);
yp=yp-mean(yp);
yar=yar-mean(yar);
yrp=yrp-mean(yrp);
yv=yv-mean(yv);
ye=ye-mean(ye);

scrsz=get(0,'ScreenSize');
pos=[0.2*scrsz(3) 0.2*scrsz(4) 0.6*scrsz(3) 0.5*scrsz(4)];

f1=figure;
%set(f,'name','Headline PCE inflation (demeaned)','numbertitle','off')
set(f1,'Position', pos);

plot(calp,yp,'-k','Linewidth',2) %inflation
hold on
plot(calp,yar+yrp,'--k','Linewidth',1) %rho
plot(calp,yv,'-k','Linewidth',1) %nu
plot(calp,ye,':k','Linewidth',1) %u

%xticks(1965:5:2000);
legend('inflation','\rho','\nu','u')
title(['Panel A. Headline inflation (demeaned), ' experiments{experiment}])

print('-bestfit',['.\figures\figure3panelA_' experiments{experiment}],'-dpdf');
%saveas(gcf, '.\figures\figure3panelA.pdf')
%[calp yp (yar+yrp) yv ye];

header={'year','inflation (all, demeaned)','rho','nu','u'};
resultsA=[header;num2cell([calp, yp, yar+yrp,yv,ye])];
xlswrite(['.\figures\figure3panelA_' experiments{experiment} '.xls'],resultsA);



% % Compute SDs of first differences %
% dyp=yp(2:end)-yp(1:end-1);
% dyv=yv(2:end)-yv(1:end-1);
% disp(['Standard Deviation of Change in PCE Inflation and V component ' num2str(std(dyp)) ', ' num2str(std(dyv))]);


% Graph a typical Series %
y=4*bpdata(:,ns_i);
x=[4*npi_ex 4*npi_fit rfac];
tmp=[calds y x];
tmp=tmp(~sum(isnan(tmp),2),:);
calp=tmp(:,1);
yp=tmp(:,2);
xp=tmp(:,3:end);
yp=yp-mean(yp);
xp=xp-mean(xp);
rp=xp(:,3:end);
v=xp(:,1);
ar=xp(:,2);
zp=yp-v-ar;
%b=zp/rp;
b=(rp'*rp)\rp'*zp;
ye=zp-rp*b;
yv=v;
yar=ar;
yrp=rp*b;
n=size(calp,1);
mis=nan(n,1);
ypma=mis;
ypma(nma+1:n)=yp(nma+1:n)-yp(1:n-nma);
yema=mis;
yema(nma+1:n)=ye(nma+1:n)-ye(1:n-nma);
yvma=mis;
yvma(nma+1:n)=yv(nma+1:n)-yv(1:n-nma);
yarma=mis;
yarma(nma+1:n)=yar(nma+1:n)-yar(1:n-nma);
yrpma=mis;
yrpma(nma+1:n)=yrp(nma+1:n)-yrp(1:n-nma);

isel = (calp > yfirst-.000001) & (calp < ylast +.0000001);
calp=calp(isel);
yp=yp(isel);
yv=yv(isel);
yar=yar(isel);
yrp=yrp(isel);
ye=ye(isel);
yp=yp-mean(yp);
yar=yar-mean(yar);
yrp=yrp-mean(yrp);
yv=yv-mean(yv);
ye=ye-mean(ye);

f2=figure;
%set(f,'name','Major household appliances inflation (demeaned)','numbertitle','off')
set(f2,'Position', pos);

plot(calp,yp,'-k','Linewidth',2)
hold on
plot(calp,yar+yrp,'--k','Linewidth',1)
plot(calp,yv,'-k','Linewidth',1)
plot(calp,ye,':k','Linewidth',1)
%xticks(1965:5:2000);
legend('inflation','\rho','\nu','u')
title(['Panel B. item ' num2str(ns_i) ' inflation (demeaned), ' experiments{experiment}])
print('-bestfit',['.\figures\figure3panelB_' experiments{experiment}],'-dpdf');
%saveas(gcf, '.\figures\figure3panelB.pdf')
%[calp yp (yar+yrp) yv ye];

header={'year',['inflation item ' num2str(ns_i) ' (demeaned)'],'rho','nu','u'};
resultsB=[header;num2cell([calp, yp, yar+yrp,yv,ye])];
xlswrite(['.\figures\figure3panelB_' experiments{experiment} '.xls'],resultsB);

% % Compute SDs of first differences %
% dyp=yp(2:end)-yp(1:end-1);
% dyv=yv(2:end)-yv(1:end-1);
% disp(['Standard Deviation of Change in PCE Inflation and V component ' num2str(std(dyp)) ', ' num2str(std(dyv))]);
% 
% % Compute Median over all series for same sample period %
% y=4*bpdata;
% yp=y(isel);
% dyp=yp(2:end,:)-yp(1:end-1,:);
% sdyp=std(dyp);
% disp(['Median across all series ' num2str(median(sdyp))]);

save(['.\figures\figure3workspace_' experiments{experiment}],'em_it_completed','npi_llf_sum','alpha','lam',...
    'rho','phi','var_eps','var_e','x3t','p3t','f','xdata','bpdata'...
    ,'nfirst','yfirst','ylast','experiment','experiments','cal_mac','d_pcep','ns_i','resultsA','resultsB');

