%SHOCKS   constructs financial shocks as in Jermann and Quadrini at an
%         annual frequency.  Run qshocks to generate quarterly forecasts.
%

%         Ellen McGrattan
%         Revised, ERM, 10/16/2016

load aggregate.csv
load aggregate_naics.csv
                   
iprint = 0;
t      = aggregate(:,1);
nt     = length(t);
x      = [aggregate(:,2:22);
          aggregate_naics(1  :45,3:23);
          aggregate_naics(46 :90,3:23);
          aggregate_naics(91 :135,3:23);
          aggregate_naics(136:180,3:23);
          aggregate_naics(181:225,3:23);
          aggregate_naics(226:270,3:23);
          aggregate_naics(271:315,3:23);
          aggregate_naics(316:360,3:23);
          aggregate_naics(361:405,3:23);
          aggregate_naics(406:450,3:23);
          aggregate_naics(451:495,3:23);
          aggregate_naics(496:540,3:23);
          aggregate_naics(541:585,3:23);
          aggregate_naics(586:630,3:23)];

label  = ['Aggregate       ';
          'Agriculture     ';
          'Mining          ';
          'Utilities       ';
          'Construction    ';
          'Manufacturing   ';
          'Wholesale trade ';
          'Retail trade    ';
          'Transportation  ';
          'Information     ';
          'FIRE            ';
          'Professional Bus';
          'Education/Health';
          'Arts/Leisure    ';
          'Other Services  '];

for i=1:15;
  xx      = x((i-1)*nt+2:i*nt,:);
  eq      = xx(:,1);
  liab    = xx(:,2);
  asset   = xx(:,3);
  netpay  = xx(:,4);
  pay     = xx(:,5);
  issue   = xx(:,6);
  div     = xx(:,7);
  erep    = xx(:,8);
  eiss    = xx(:,9);
  int     = xx(:,10);
  ldrep   = xx(:,11);
  ldiss   = xx(:,12);
  sdnet   = xx(:,13);
  sale    = xx(:,14);
  capx    = xx(:,15);
  aqc     = xx(:,16);
  sppe    = xx(:,17);
  oancf   = xx(:,18);
  fopt    = xx(:,19);
  cpi     = xx(:,20);
  infl    = xx(:,21);

  eqpay   = div+erep-eiss;
  dbpay   = int+ldrep-ldiss+sdnet;
  deld    = ldrep-ldiss+sdnet;

  Eq(:,i)     = eq;
  Liab(:,i)   = liab;
  Capx(:,i)   = capx;
  Sale(:,i)   = sale;
  Fopt(:,i)   = fopt;
  Oancf(:,i)  = oancf;
  Aqc(:,i)    = aqc;
  Sppe(:,i)   = sppe;
  X(:,i)      = capx+aqc-sppe;
  K0          = eq(1);
  for s=1:nt-1;
    Keoy(s,i) = (1-.032)*K0+X(s,i);
    K0        = Keoy(s,i);
  end;
  Phi(:,i)    = sale./(Keoy(:,i)-liab)/4;
end;   
if (iprint==1);
  for i=1:15;
    figure(1)
    plot(t(2:nt),Keoy(:,i)./Sale(:,i),t(2:nt),Liab(:,i)./Sale(:,i))
    title(label(i,:))
    legend('Capital/Sales','Liabilities/Sales')
    axis([1972,2015,0,2.5])
    figure(2)
    plot(t(2:nt),Phi(:,i))
    title(label(i,:))
    pause(1)
  end;
end;

shks  = Phi(:,[3,6,9,14]);
figure,plot(t(15:nt),shks(14:nt-1,:)),
title('Financial Shocks')
legend('Mining','Manufacturing','Transportation','Arts and Leisure')

tab = [t(2:nt),shks];

