clear all
close all
clc
warning off;

%% EDITED PROGRAM FROM RAMEY-ZUBAIRY
%% DOUBLE PERCENT SIGNS ARE ANNOTATIONS FOR THOSE WHO ARE NEW TO MATLAB

%% The ; tells it not to display

%% Read the data from the Excel file.
%% The "1" below refers to the worksheet

vdata = xlsread('RZ_Data_2013Sept30.xlsx',1); % US data

%% time is a new variable,
%% the 1 is the first column of the vdata

time=vdata(:,1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%IMPORTANT INPUT (sample window)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% find will look in vector for time==1890

t1=find(time==1890);
t2=find(time==2010.75);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%IMPORTANT INPUT (threshold for the unemployment rate)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
threshold=6.5;

%% The following selects particular rows of whole vdata matrix
%% The : means all of the columns
%% We could have said vdatasub = vdata(t1:t2,:)


vdata=vdata(t1:t2,:); 

%% q is the quarter, now choose the first column but all the rows
%% First argument is row, second is column
%% Each variable is a vector we are extracting

q=vdata(:,1);
ngov=vdata(:,2);
ngdp=vdata(:,3);
pgdp=vdata(:,4);
totpop=vdata(:,5);
recession=vdata(:,6);
unemp=vdata(:,7);
news=vdata(:,8);
realgdp=vdata(:,9);
ntax=vdata(:,10);
rgdp_pot=vdata(:,11);


%Creating variables of interest

%% Since news can take negative values, we can't take logs
%% Thus, scale it by lagged nominal GDP

%% Start with t=2 because know we need to divide by lag of GDP

%% For changing one value, we could have newsy(1,1) = 0;

for t=2:size(news,1)
newsy(t,1)= news(t,1)./ngdp(t-1,1);
end


%% Use . when you need to multiply, divide, etc. element by element;

rgdp=(realgdp./totpop);
rgov=(ngov./totpop./pgdp);

lrgdp=log(realgdp./totpop);
lrgov=log(ngov./totpop./pgdp);


%% Find unemployment where >= threshold value, say element 3, 6, 8, 9

rrr=find(unemp>=threshold); 
rrrn=rrr+1; %Since unemployment above threshold in the previous period is the criterion

%% length is how many elements;

if rrrn(end)==length(unemp)+1
    rrrn=rrrn(1:end-1);
else
    rrrn=rrrn;
end

%% fu is a vector of zeroes that is the size of unemployment
fu=zeros(size(unemp));

%% replacing the 0's with 1's when U(t-1) >= threshold

fu(rrrn)=1; % This is the dummy variable that defines the two states


nlags=5; %Number of lags we are interested in +1

%% Define lags correctly

y=lrgdp(nlags:end);
fu=fu(nlags:end);

Ly=lrgdp(nlags-1:end-1);
L2y=lrgdp(nlags-2:end-2);
L3y=lrgdp(nlags-3:end-3);
L4y=lrgdp(nlags-4:end-4);

Lg=lrgov(nlags-1:end-1);
L2g=lrgov(nlags-2:end-2);
L3g=lrgov(nlags-3:end-3);
L4g=lrgov(nlags-4:end-4);

%% vector of ones for constant term

constant= ones(length(y),1);

%% start at 1, increment of 1, to length of y

t=1:1:length(y);
tsq=t.^2;
tcu=t.^3;
tqu=t.^4;

%Creating part of the new variables that enter the RHS, which are in
%levels, not logs.
ynew=rgdp(nlags:end);
Lynew=rgdp(nlags-1:end-1);

gnew=rgov(nlags:end);
Lgnew=rgov(nlags-1:end-1);

ramey=newsy(nlags:end); % shock of interest

%% horizon for impulse response

hor=20;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%LINEAR
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% t' tsq' etc. is the transpose - was created as row, want to turn into column;

x=[constant, t', tsq', tcu', tqu', ramey, Lg, L2g, L3g, L4g, Ly, L2y, L3y, L4y];

rpos=6; %position of the shock

%% x is the whole matrix with regressors, r is the # of rows, nn is the number of columns;

[r,nn]=size(x);

%% Loop of regressions
%% if NA, Matlab can't deal with it like Stata does

%% i, 1 because of Jorda type stuff - last part ensures that vector is of the same length.

for i=1:hor
    yy=(ynew(i:end)- Lynew(1:end-i+1))./Lynew(1:end-i+1); %Y variable entering in the LHS

%% Calls program nwest, first argument is yy second is matrix x, need to make sure dimensions are right. first part is rows, after "," is columns
%%  i adjusts for horizon.
%% Results is specific to nwest program.  Output is called results.

    results=nwest(yy, x(1:end-i+1,:),i);

%% how coefficients are stored.
    b=results.beta;
    bint(:,1)=b-(results.se*1.96);
    bint(:,2)=b+(results.se*1.96);

%% how standard errors are stored
    se(:,i)=results.se';

%% reglin stores everything
    reglin(:,i)=b;
    confidencelin(:,:,i)=bint;
end

%% confidence bands - 3 dimensional thing

confidencelin=reshape(confidencelin,2*nn,hor);

%% extracting the coefficient to ramey variable

linGDP=[reglin(rpos,:);se(rpos,:)];

clear reg confidence bint se

%% still in linear, now running with gov instead y - just redefined yy
y=lrgov(nlags:end);
for i=1:hor
    yy=(gnew(i:end)- Lgnew(1:end-i+1))./Lynew(1:end-i+1); %G variable entering in the LHS
    results=nwest(yy, x(1:end-i+1,:),i);
    b=results.beta;
    bint(:,1)=b-(results.se*1.96);
    bint(:,2)=b+(results.se*1.96);
    reg(:,i)=b;
    se(:,i)=results.se';
    confidence(:,:,i)=bint;
end
confidence=reshape(confidence,2*nn,hor);
linGOV=[reg(rpos,:);se(rpos,:)];

%% vector of zeroes for figures - need 0 axis

zz=zeros(1,hor);

%% Figure is one row, 2 columns
%% figure(3) is what the figure will be called

figure(3)

%% can do simple too:  plot time,unemp

%% need two different subplots in figure, 1 row, 2 columns, 1st of the subplot figures
subplot(1,2,1)

%% plot is command for figure, x axis is 1 to horizon with increment of 1, coefficients of rpos

plot(1:1:hor, reg(rpos,:))

% hold back and then also plot something else on same figure

hold on

%% this plots the zeroes as well on the same figure

plot(1:1:hor, zz, 'k-')
hold on
grpyat=[(1:1:hor)', confidence(rpos,:)'; (hor:-1:1)' confidence(rpos+nn,hor:-1:1)'];
patch(grpyat(:,1), grpyat(:,2), [0.7 0.7 0.7],'edgecolor', [0.7 0.7 0.7]);
plot(1:1:hor, zz, 'k-')
plot(1:1:hor, reg(rpos,:), 'k','LineWidth', 1.5)
title('Government spending');
xlabel('quarter')
axis tight
subplot(1,2,2)
plot(1:1:hor, reglin(rpos,:))
hold on
plot(1:1:hor, zz, 'k-')
hold on
grpyat=[(1:1:hor)', confidencelin(rpos,:)'; (hor:-1:1)' confidencelin(rpos+nn,hor:-1:1)'];
patch(grpyat(:,1), grpyat(:,2), [0.7 0.7 0.7],'edgecolor', [0.7 0.7 0.7]);%[0.65 0.65 0.65]);
plot(1:1:hor, zz, 'k-')
plot(1:1:hor, reglin(rpos,:), 'k','LineWidth', 1.5)
title('GDP');
xlabel('quarter')
axis tight



% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %NON-LINEAR
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear reg b bint se
clear confidence

x=[constant, t', tsq', tcu', tqu', fu, fu.*ramey, (1-fu).*ramey, fu.*Lg, fu.*L2g, fu.*L3g, fu.*L4g, fu.*Ly, fu.*L2y, fu.*L3y, fu.*L4y, ...
                                                (1-fu).*Lg, (1-fu).*L2g, (1-fu).*L3g, (1-fu).*L4g, (1-fu).*Ly, (1-fu).*L2y, (1-fu).*L3y, (1-fu).*L4y];
[r,nn]=size(x);
rpos=7;%position of the shock

for i=1:hor
    yy=(ynew(i:end)- Lynew(1:end-i+1))./Lynew(1:end-i+1); %Y variable entering in the LHS
    results=nwest(yy, x(1:end-i+1,:),i);
    b=results.beta;
    bint(:,1)=b-(results.se*1.96);
    bint(:,2)=b+(results.se*1.96);
    se(:,i)=results.se';
    regy(:,i)=b;
    confidencey(:,:,i)=bint;
end
confidencey=reshape(confidencey,2*nn,hor);
stateaGDP=[regy(rpos,:);se(rpos,:)];
statebGDP=[regy(rpos+1,:);se(rpos+1,:)];


clear reg b bint se
clear confidence


x=[constant, t', tsq', tcu', tqu', fu, fu.*ramey, (1-fu).*ramey, fu.*Lg, fu.*L2g, fu.*L3g, fu.*L4g, fu.*Ly, fu.*L2y, fu.*L3y, fu.*L4y, ...
                                                (1-fu).*Lg, (1-fu).*L2g, (1-fu).*L3g, (1-fu).*L4g, (1-fu).*Ly, (1-fu).*L2y, (1-fu).*L3y, (1-fu).*L4y];
[r,nn]=size(x);


for i=1:hor
    yy=(gnew(i:end)- Lgnew(1:end-i+1))./Lynew(1:end-i+1); %G variable entering in the LHS
    results=nwest(yy, x(1:end-i+1,:),i);
    b=results.beta;
    bint(:,1)=b-(results.se*1.96);
    bint(:,2)=b+(results.se*1.96);
    se(:,i)=results.se';
    reg(:,i)=b;
    confidence(:,:,i)=bint;
end
confidence=reshape(confidence,2*nn,hor);
stateaGOV=[reg(rpos,:);se(rpos,:)];
statebGOV=[reg(rpos+1,:);se(rpos+1,:)];


figure(4)
subplot(1,2,1)
plot(1:1:hor, reg(rpos,:),'k')
hold on
grpyat=[(1:1:hor)', confidence(rpos,:)'; (hor:-1:1)' confidence(rpos+nn,hor:-1:1)'];
patch(grpyat(:,1), grpyat(:,2), [0.7 0.7 0.7],'edgecolor', [0.7 0.7 0.7]);
plot(1:1:hor, reg(rpos,:), 'k','LineWidth', 1.5)
hold on
plot(1:1:hor, reg(rpos+1,:), 'r-o', 1:1:hor, confidence(rpos+1,:), 'r--', 1:1:hor, confidence(rpos+1+nn,:), 'r--', 'LineWidth', 1.5);
axis tight
title('Government Spending')
xlabel('quarter')
subplot(1,2,2)
plot(1:1:hor, regy(rpos,:),'k')
hold on
grpyat=[(1:1:hor)', confidencey(rpos,:)'; (hor:-1:1)' confidencey(rpos+nn,hor:-1:1)'];
patch(grpyat(:,1), grpyat(:,2), [0.7 0.7 0.7],'edgecolor', [0.7 0.7 0.7]);
plot(1:1:hor, regy(rpos,:), 'k','LineWidth', 1.5)
hold on
plot(1:1:hor, regy(rpos+1,:), 'r-o', 1:1:hor, confidencey(rpos+1,:), 'r--', 1:1:hor, confidencey(rpos+1+nn,:), 'r--', 'LineWidth', 1.5);
axis tight
title('GDP')
xlabel('quarter')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Multipliers:

topasteinExcel=[linGDP;linGOV;stateaGDP;statebGDP;stateaGOV;statebGOV];
disp ('Output multipliers for the US')
disp ( 'The multipliers by row are: 2 year integral multiplier, 4 year integral multiplier and peak multiplier')
disp('LINEAR ')
[sum(linGDP(1,1:8))./sum(linGOV(1,1:8)), sum(linGDP(1,1:16))./sum(linGOV(1,1:16)), max(linGDP(1,1:16))./max(linGOV(1,1:16))]
disp('HIGH UNEMPLOYMENT state')
[sum(stateaGDP(1,1:8))./sum(stateaGOV(1,1:8)), sum(stateaGDP(1,1:16))./sum(stateaGOV(1,1:16)), max(stateaGDP(1,1:16))./max(stateaGOV(1,1:16))]
disp('LOW UNEMPLOYMENT state')
[sum(statebGDP(1,1:8))./sum(statebGOV(1,1:8)), sum(statebGDP(1,1:16))./sum(statebGOV(1,1:16)), max(statebGDP(1,1:16))./max(statebGOV(1,1:16))]