%RESULTS  prints out results in:
% 
%           Intangible Capital and Measured Productivity
%           Minneapolis Fed Staff Report 
%


%----------------------------------------------------
% Variance Decomposition with Ergodic Distribution
%----------------------------------------------------
load ../MLEdata/ldhpctot.dat
load ../MLEdata/ldgypc.dat
load ../MLEdata/ldhpc21.dat
load ../MLEparam/fixedpar.dat
load ../MLEest/mleest.out

global ZVAR
ZVAR  = [ldgypc(149:268,[1:15]),ldhpc21(149:268,[20,17,19])];
[L,innov,Xt,G,H,J,J2,XT] = mleq(mleest,fixedpar);
Sig         = H*H';
nz          = [31:60,78,82];
[tab,tab1]  = vardec(G(1:82,1:82),Sig(1:82,1:82)+1e-10*eye(82),nz,J(:,1:82));
for i=1:15;
  tab2(i,2) = tab1(i,i)+tab1(i,i+15);
  tab2(i,3) = sum(tab1(i,1:30))-tab2(i,2);
  tab2(i,1) = tab2(i,2)+tab2(i,3);
end;
tab2(16,2)  = tab1(16,5)+tab1(16,20);
tab2(16,3)  = sum(tab1(16,1:30))-tab2(16,2);
tab2(16,1)  = tab2(16,2)+tab2(16,3);
tab2(17,2)  = tab1(17,9)+tab1(17,24);
tab2(17,3)  = sum(tab1(17,1:30))-tab2(17,2);
tab2(17,1)  = tab2(17,2)+tab2(17,3);
tab2(18,2)  = tab1(18,11)+tab1(18,26);
tab2(18,3)  = sum(tab1(18,1:30))-tab2(18,2);
tab2(18,1)  = tab2(18,2)+tab2(18,3);
for i=1:18;
  tab2(i,4) = tab1(i,32);
  tab2(i,5) = tab1(i,31);
end;
tab2        = tab2*100;

disp(' ')
disp(' ')
disp('  Variance Decomposition using the Ergodic Distribution')
disp('-----------------------------------------------------------')
disp('                             TFP shocks             | Labor')
disp('Observable:        Sectoral   Own    Other  Common  | Wedge')
disp('-----------------------------------------------------------')
disp('Gross output in:')
disp(sprintf('  Agriculture      %7.1f %7.1f %7.1f %7.1f %7.1f',tab2(1,:)))
disp(sprintf('  Mining           %7.1f %7.1f %7.1f %7.1f %7.1f',tab2(2,:)))
disp(sprintf('  Utilities        %7.1f %7.1f %7.1f %7.1f %7.1f',tab2(3,:)))
disp(sprintf('  Construction     %7.1f %7.1f %7.1f %7.1f %7.1f',tab2(4,:)))
disp(sprintf('  Manufacturing    %7.1f %7.1f %7.1f %7.1f %7.1f',tab2(5,:)))
disp(sprintf('  Wholesale trade  %7.1f %7.1f %7.1f %7.1f %7.1f',tab2(6,:)))
disp(sprintf('  Retail trade     %7.1f %7.1f %7.1f %7.1f %7.1f',tab2(7,:)))
disp(sprintf('  Transportation   %7.1f %7.1f %7.1f %7.1f %7.1f',tab2(8,:)))
disp(sprintf('  Information      %7.1f %7.1f %7.1f %7.1f %7.1f',tab2(9,:)))
disp(sprintf('  FIRE             %7.1f %7.1f %7.1f %7.1f %7.1f',tab2(10,:)))
disp(sprintf('  Prof & business  %7.1f %7.1f %7.1f %7.1f %7.1f',tab2(11,:)))
disp(sprintf('  Education etc    %7.1f %7.1f %7.1f %7.1f %7.1f',tab2(12,:)))
disp(sprintf('  Leisure etc.     %7.1f %7.1f %7.1f %7.1f %7.1f',tab2(13,:)))
disp(sprintf('  Other services   %7.1f %7.1f %7.1f %7.1f %7.1f',tab2(14,:)))
disp('Hours in:')
disp(sprintf('  Comp & Elec      %7.1f %7.1f %7.1f %7.1f %7.1f',tab2(16,:)))
disp(sprintf('  BC & Telecomm    %7.1f %7.1f %7.1f %7.1f %7.1f',tab2(17,:)))
disp(sprintf('  Advertising      %7.1f %7.1f %7.1f %7.1f %7.1f',tab2(18,:)))
disp('-----------------------------------------------------------')
disp(' ')
disp(' ')


%----------------------------------------------------
% Variance Decomposition with In-Sample Projections
%----------------------------------------------------
hss         = .262;
hact        = ldhpctot(149:267);
hact        = exp(hact-mean(hact))*hss;
yact        = ldgypc(149:267,22)-mean(ldgypc(149:267,22));
ypred       = XT*J2(2,:)'-mean(XT*J2(2,:)');

hpred       = XT*J2(1,:)';
hpred       = exp(hpred-mean(hpred))*hss;

ii          = [1:30,61:83];
hpredc      = XT(:,ii)*J2(1,ii)';
hpredc      = exp(hpredc-mean(hpredc))*hss;
ypredc      = XT(:,ii)*J2(2,ii)'-mean(XT(:,ii)*J2(2,ii)');

ii          = [1:81,83];
hpreds      = XT(:,ii)*J2(1,ii)';
hpreds      = exp(hpreds-mean(hpreds))*hss;
ypreds      = XT(:,ii)*J2(2,ii)'-mean(XT(:,ii)*J2(2,ii)');

t           = ptime(1985.2,119,4);
figure(1)
plot(t,hact,t,hpred)
title('Per Capita Hours of Work, 1985:1-2014:4')
legend('Actual','Predicted')

figure(2)
plot(t,yact,t,ypred)
title('Logged and Detrended Gross Output, 1985:1-2014:4')
legend('BEA','Predicted')

figure(3)
plot(t,hpred,t,hpredc,t,hpreds)
title('Predicted Per Capita Hours of Work, 1985:1-2014:4')
legend('All Shocks','Only Common','Only Sectoral','Location','SouthEast')

figure(4)
plot(t,ypred,t,ypredc,t,ypreds)
title('Predicted Logged and Detrended Gross Output, 1985:1-2014:4')
legend('All Shocks','Only Common','Only Sectoral','Location','SouthWest')

hfigs       = [t,hact,hpred,hpredc,hpreds];
yfigs       = [t,yact,ypred,ypredc,ypreds];

%----------------------------------------------------
% Variance decomposition for 1990s tech boom
%----------------------------------------------------

t1          = find(t==1991.875);
t2          = find(t==2000.625);
for i=1:16;
   ii       = 30+i;
   if (i==5);
     ii     = [35,50];
   end;
   if (i==9);
     ii     = [39,54];
   end;
   if (i==11);
     ii     = [41,56];
   end;
   if (i==16);
     ii     = 82;
   end;
   dh(i,1)  = (XT(t2,ii)*J2(1,ii)'-XT(t1,ii)*J2(1,ii)')*100;
   dy(i,1)  = (XT(t2,ii)*J2(2,ii)'-XT(t1,ii)*J2(2,ii)')*100;
end;

%----------------------------------------------------
% Variance decomposition for Great Recession
%----------------------------------------------------

t1          = find(t==2007.875);
t2          = find(t==2009.625);
for i=1:16;
   ii       = 30+i;
   if (i==5);
     ii     = [35,50];
   end;
   if (i==9);
     ii     = [39,54];
   end;
   if (i==11);
     ii     = [41,56];
   end;
   if (i==16);
     ii     = 82;
   end;
   dh(i,2)  = (XT(t2,ii)*J2(1,ii)'-XT(t1,ii)*J2(1,ii)')*100;
   dy(i,2)  = (XT(t2,ii)*J2(2,ii)'-XT(t1,ii)*J2(2,ii)')*100;
end;
figure(5),bar(dy)
xlabel('15 Major Sectors and Common Shock')
title('Decomposition of Changes in Gross Output in Two Periods')
legend('Tech Boom','Great Recession','Location','Northwest')

disp('  Decomposition of Changes in Gross Output, Two Periods ')
disp('-----------------------------------------------------------')
disp('                           Technology      Great ')
disp('                              Boom        Recession')
disp('-----------------------------------------------------------')
disp(sprintf('Common shock         %13.1f %13.1f ',dy(16,:)))
disp('Sectoral shocks: ')          
disp(sprintf('  Agriculture        %13.1f %13.1f ',dy(1,:)))
disp(sprintf('  Mining             %13.1f %13.1f ',dy(2,:)))
disp(sprintf('  Utilities          %13.1f %13.1f ',dy(3,:)))
disp(sprintf('  Construction       %13.1f %13.1f ',dy(4,:)))
disp(sprintf('  Manufacturing      %13.1f %13.1f ',dy(5,:)))
disp(sprintf('  Wholesale trade    %13.1f %13.1f ',dy(6,:)))
disp(sprintf('  Retail trade       %13.1f %13.1f ',dy(7,:)))
disp(sprintf('  Transportation     %13.1f %13.1f ',dy(8,:)))
disp(sprintf('  Information        %13.1f %13.1f ',dy(9,:)))
disp(sprintf('  FIRE               %13.1f %13.1f ',dy(10,:)))
disp(sprintf('  Prof & business    %13.1f %13.1f ',dy(11,:)))
disp(sprintf('  Education etc      %13.1f %13.1f ',dy(12,:)))
disp(sprintf('  Leisure etc.       %13.1f %13.1f ',dy(13,:)))
disp(sprintf('  Other services     %13.1f %13.1f ',dy(14,:)))
disp('-----------------------------------------------------------')
disp(sprintf('Total change         %13.1f %13.1f ',[sum(dy(:,1)),sum(dy(:,2))]))
disp('-----------------------------------------------------------')
disp(' ')
disp(' ')




%----------------------------------------------------
% Compare US Solow Residual and Model's Common TFP
%----------------------------------------------------

figure(6)
load ../TFPdata/lzt.dat
lzt  = .67*lzt(2:120);
tem1 = (lzt-lzt(91))/std(lzt);
tem2 = (XT(:,82)-XT(91,82))/std(XT(:,82));
plot(t,[tem1,tem2])
title('Solow Residual and Common TFP, 1985:1-2014:4')
legend('Solow Residual','Common TFP')
c1 = corrcoef([tem1,tem2]);
c2 = corrcoef([tem1-hptrend(tem1,1600),tem2-hptrend(tem2,1600)]);
text(1988,-2.0,['Raw correlation = ',num2str(c1(1,2)*100),'%'])
text(1988,-2.3,['HP correlation  = ',num2str(c2(1,2)*100),'%'])

figure(7)
plot(t,[tem1,tem2])
title('Solow Residual and Common TFP, 1985:1-2014:4')
legend('Solow Residual','Common TFP')
axis([2007.875,2015,-3,1])

%----------------------------------------------------
% Compare Fernald's TFP and Model's Common TFP
%----------------------------------------------------

figure(8)
load ../TFPdata/rawdata/Fernald/qtrly.dat
lfer   = cumsum(qtrly(153:271,15))/400;
Xfer   = [ones(119,1),[1:119]'];
beta   = ols(Xfer,lfer);
lfer   = lfer-Xfer*beta;
tem3   = (lfer-lfer(91))/std(lfer);
plot(t,[tem3,tem2])
title('Fernald TFP and Common TFP, 1985:1-2014:4')
legend('Fernald','Common TFP')
c1 = corrcoef([tem3,tem2]);
c2 = corrcoef([tem3-hptrend(tem3,1600),tem2-hptrend(tem2,1600)]);
text(1988,-2.0,['Raw correlation = ',num2str(c1(1,2)*100),'%'])
text(1988,-2.3,['HP correlation  = ',num2str(c2(1,2)*100),'%'])

figure(9)
plot(t,[tem3,tem2])
title('Fernald TFP and Common TFP, 1985:1-2014:4')
legend('Fernald','Common TFP')
axis([2007.875,2015,-3,1])

%----------------------------------------------------
% Sectoral Implications that serve as external checks
%----------------------------------------------------

dgytot  = ldgypc(149:267,22)-hptrend(ldgypc(149:267,22),1600);
dgytotp = ypred-hptrend(ypred,1600);
dgyact  = ldgypc(149:267,1:15)-hptrend(ldgypc(149:267,1:15),1600);
dgypred = XT*J2(3:17,:)'-hptrend(XT*J2(3:17,:)',1600);

dhact   = ldhpc21(149:267,1:15)-hptrend(ldhpc21(149:267,1:15),1600);
dhpred  = XT*J2(18:32,:)'-hptrend(XT*J2(18:32,:)',1600);
dh1pred = XT*J2(33:47,:)'-hptrend(XT*J2(33:47,:)',1600);
dh2pred = XT*J2(48:50,:)'-hptrend(XT*J2(48:50,:)',1600);

dpract   = ldgypc(149:267,1:15)-ldhpc21(149:267,1:15);
dpract   = dpract-hptrend(dpract,1600);
dprpred  = XT*J2(3:17,:)'-XT*J2(18:32,:)';
dprpred  = dprpred-hptrend(dprpred,1600);


ii      = 1:15;
[c,s,l] = pca(dhact(:,ii)./(ones(119,1)*std(dhact(:,ii))));
fl(:,1) = c(:,1)*100;
apca    = l(1)/sum(l)*100;
[c,s,l] = pca(dhpred(:,ii)./(ones(119,1)*std(dhpred(:,ii))));
fl(:,2) = c(:,1)*100;
ppca    = l(1)/sum(l)*100;

disp('     Principal Component Analysis of Sectoral Hours        ')
disp('-----------------------------------------------------------')
disp('                               Actual      Predicted')
disp('-----------------------------------------------------------')
disp('Factor loadings: ')
disp(sprintf('  Agriculture        %13.0f %13.0f',fl(1,:)))
disp(sprintf('  Mining             %13.0f %13.0f',fl(2,:)))
disp(sprintf('  Utilities          %13.0f %13.0f',fl(3,:)))
disp(sprintf('  Construction       %13.0f %13.0f',fl(4,:)))
disp(sprintf('  Manufacturing      %13.0f %13.0f',fl(5,:)))
disp(sprintf('  Wholesale trade    %13.0f %13.0f',fl(6,:)))
disp(sprintf('  Retail trade       %13.0f %13.0f',fl(7,:)))
disp(sprintf('  Transportation     %13.0f %13.0f',fl(8,:)))
disp(sprintf('  Information        %13.0f %13.0f',fl(9,:)))
disp(sprintf('  FIRE               %13.0f %13.0f',fl(10,:)))
disp(sprintf('  Prof & business    %13.0f %13.0f',fl(11,:)))
disp(sprintf('  Education etc      %13.0f %13.0f',fl(12,:)))
disp(sprintf('  Leisure etc.       %13.0f %13.0f',fl(13,:)))
disp(sprintf('  Other services     %13.0f %13.0f',fl(14,:)))
disp('-----------------------------------------------------------')
disp(sprintf('PC variance          %13.0f %13.0f',[apca,ppca]))
disp('-----------------------------------------------------------')


corrs   = zeros(14,2);
for i=1:14;
  tem        = corrcoef([dpract(:,i),dprpred(:,i),dhact(:,i),dhpred(:,i)]);
  corrs(i,:) = [tem(1,2),tem(3,4)]*100;
end;
wgt     = [0.0063;0.0097;0.0126;0.0488;0.1436;0.0746;0.0784;
           0.0396;0.0326;0.1159;0.1946;0.1458;0.0541;0.0433]; %=share of hours

disp(' ')
disp(' ')
disp('  Correlations between Predicted and Actual Sectoral Data')
disp('-----------------------------------------------------------')
disp('                            Productivity    Hours ')
disp('-----------------------------------------------------------')
disp(sprintf('  Agriculture        %13.0f %13.0f',corrs(1,:)))
disp(sprintf('  Mining             %13.0f %13.0f',corrs(2,:)))
disp(sprintf('  Utilities          %13.0f %13.0f',corrs(3,:)))
disp(sprintf('  Construction       %13.0f %13.0f',corrs(4,:)))
disp(sprintf('  Manufacturing      %13.0f %13.0f',corrs(5,:)))
disp(sprintf('  Wholesale trade    %13.0f %13.0f',corrs(6,:)))
disp(sprintf('  Retail trade       %13.0f %13.0f',corrs(7,:)))
disp(sprintf('  Transportation     %13.0f %13.0f',corrs(8,:)))
disp(sprintf('  Information        %13.0f %13.0f',corrs(9,:)))
disp(sprintf('  FIRE               %13.0f %13.0f',corrs(10,:)))
disp(sprintf('  Prof & business    %13.0f %13.0f',corrs(11,:)))
disp(sprintf('  Education etc      %13.0f %13.0f',corrs(12,:)))
disp(sprintf('  Leisure etc.       %13.0f %13.0f',corrs(13,:)))
disp(sprintf('  Other services     %13.0f %13.0f',corrs(14,:)))
disp('-----------------------------------------------------------')
disp(sprintf('  Weighted Avg       %13.0f %13.0f', ...
                [sum(corrs(1:14,1).*wgt),sum(corrs(1:14,2).*wgt)]))
disp('-----------------------------------------------------------')
disp(' ')
disp(' ')


%----------------------------------------------------
% Implications for latent TFPs and Intangibles
%----------------------------------------------------

tfps    = XT(:,[31:35,50,36:39,54,40:41,56,42:44,82]);
xtu     = XT*J2(51:65,:)';
xiu     = XT*J2(66:80,:)';
xtc     = XT*J2(81:95,:)';
xic     = XT*J2(96:110,:)';
xttot   = XT*J2(111,:)';
xitot   = XT*J2(112,:)';

labels  = ['TFP, Agriculture       '
           'TFP, Mining            '
           'TFP, Utilities         '
           'TFP, Construction      '
           'TFP, Manufacturing     '
           'TFP, Manuf, Comp&Elec  '
           'TFP, Wholesale trade   '
           'TFP, Retail trade      '
           'TFP, Transportation    '
           'TFP, Information       '
           'TFP, Info, BC&Telecomm '
           'TFP, FIRE              '
           'TFP, Prof & Business   '
           'TFP, P&B, Advertising  '
           'TFP, Education etc     '
           'TFP, Leisure etc       '
           'TFP, Other services    '
           'TFP, Common            '
           'XIC, Manufacturing     '
           'XIC, Information       '
           'XIC, Prof & Business   '];

tab     = kpstats([tfps,xic(:,[5,9,11])],XT*J2(2,:)',1,labels);

figure(10)
load ../TFPdata/tanginv.dat
plot(t,xitot/xitot(91)*100,t,xttot/xttot(91)*100,tanginv(:,1),tanginv(:,2));
title('Predicted and Actual Investments, 1985:1-2014:4')
legend('Predicted Intangible','Predicted Tangible','US Tangible')
axis([2007.875,2015,50,110])
