
% This m-file provides a template for estimating "spatial impulse response
% functions," as discussed in Brady (2009). I say template since the spatial
% IRFs are estimated by cobbling together different methods--estimating Jorda's
% (2005) linear projections from a spatial autoregressive model.  

% The spatial code used in this paper is from James Lesage (see his website
% www.spatial-econometrics.com for the spatial code referenced in this
% file). 

% In this template I attempt to provide a guide to putting the methods
% together.  Naturally, one will have to have to do their own leg work in
% getting their data ready--however, where applicable I discuss my own
% application to make it easier to "cut and paste" your own data.

% If you have any questions, or find errors, please email me, Ryan Brady, at
% rbrady@usna.edu.  

%References:
% Brady, Ryan, R., (2009), "Measuring the diffusion of housing prices
% across space and over time, Journal of Applied Econometrics, 
% DOI:10.1002/jae.1118
% Jorda, Oscar, (2005), "Estimation and inferences of impulse responses by
% local projections. American Economic Review 95: 61?82.


%1) Getting the Data prepared: Note that I prepared much of my data in
%Excel and Stata--then I imported the data into Matlab. For those who are
%more skilled in Matlab than I was when I estimated the spatial IRFs, you
%probably won't need much of the crude coding in this section.

% My commands: 
clear
load MSA %my prepared data loaded from the Matlab workspace (previously saved as a .mat file)

T = 179;  %Time in months for each region-1
nc = 26; %The number of regions

reg = data(:,1); %County numbers
dt = data(:,4); %date
y = data(:,16); %housing prices (natural log)
y1 = data(:,17); %housing prices lag
x1 = data(:,8);  %unemployment
x2 = data(:,21);  %construction (natural log)
x3 = data(:,10);  %realrate
x4 = data(:,22);  %population (natural log)
x5 = data(:,23); %industrial production (natural log)
x6 = data(:,12:14);  %season dummies
x7 = data(:,44:68);  %county dummies (46 for 47 counties)

%lags of the regressors
y2 = data(:,18); y3 = data(:,19); y4 = data(:,20); %more housing lags
x1_1 = data(:,28); x1_2 = data(:,29); x1_3 = data(:,30); x1_4 = data(:,31); %unemployment lags
x2_1 = data(:,24); x2_2 = data(:,25); x2_3 = data(:,26); x2_4 = data(:,27); %construction lags
x3_1 = data(:,40); x3_2 = data(:,41); x3_3 = data(:,42); x3_4 = data(:,43); %realrate lags
x4_1 = data(:,32); x4_2 = data(:,33); x4_3 = data(:,34); x4_4 = data(:,35); %population lags
x5_1 = data(:,36); x5_2 = data(:,37); x5_3 = data(:,38); x5_4 = data(:,39); %IP lags
yleads = data(:,69:128);  %Housing price leads +1 through +60 (the third column in leads is the current hprice)

% 2) Here I define the initial weighting matrix.  James Lesage does provide
% code for automating this part, but I chose to define the matrix by
% "hand."

% weighting matrix is from excel file

% In the next few lines I sort the data
lags = [y2 y3 y4 x1_1 x1_2 x1_3 x1_4 x2_1 x2_2 x2_3 x2_4 x3_1 x3_2 x3_3 x3_4 x4_1 x4_2 x4_3 x4_4 x5_1 x5_2 x5_3 x5_4];
allp = [dt reg y y1 x1 x2 x3 x4 x5 x6 x7 yleads lags];
[n,m] = size(allp);

% Here I rid the data matrix "all" of the observations lost due to the lags
% used in estimation (i.e., if you decide on an AR(1), or AR(2) . . . )
for k = 1:n;
     %if k == 2671, break, end  %This should be plus 1 of the final sample (after taking lags)
 if k == 4655, break, end
 if allp(k,1) == 1; %Since I sorted the data in Excel, I had some zeros (after taking lags), which I changed to .0012 to import into Matlab
    allp(k,:) = [];        %This line gets rid of rows with any cell equal to 0.0012.
 end;
end;


allp = sortrows(allp,1); %This sorts the matrix according to data ("dt"), this is required by Lesage's code used below.
[n,m] = size(allp); 
y =  (allp(:,3))*100; %Here I just redefine the variables in the "all" matrix
y1 = (allp(:,4))*100;
x1 = allp(:,5);
x2 = (allp(:,6))*100;
x3 = allp(:,7);
x4 = (allp(:,8))*100;
x5 = (allp(:,9))*100;
x6 = allp(:,10:12);
x7 = allp(:,13:37);
dt = allp(:,1);
yleads = (allp(:,38:97))*100; %36 Leads
lags = (allp(:,98:m))*100;
W = x; %Rename the weighting matrix consistent with Lesage's code.

% Lesage's code: "sparse" "normw" and "slag" are from Lesage's spatial toolbox (see his webiste).  
WW = sparse(W);
W = normw(WW);
W = slag(W,2);

% This code creates the "spatial lag" variable denoted in the literature as Wy.
% I also estimate the "lags" for the regressors for use in IV.
% These commands are based on Paul Elhorst's "sar_panel.m" file
% provided in Lesage's toolkit.

for t=1:T
    t1=1+(t-1)*nc;t2=t*nc;
    Wy([t1:t2],1)= sparse(W)*y([t1:t2],1);
end
%Instruments
for t=1:T
    t1=1+(t-1)*nc;t2=t*nc;
    Wx1([t1:t2],1)= sparse(W)*x1([t1:t2],1);
end
for t=1:T
    t1=1+(t-1)*nc;t2=t*nc;
    Wx2([t1:t2],1)= sparse(W)*x2([t1:t2],1);
end
for t=1:T
    t1=1+(t-1)*nc;t2=t*nc;
    Wx4([t1:t2],1)= sparse(W)*x4([t1:t2],1);
end
for t=1:T
    t1=1+(t-1)*nc;t2=t*nc;
    Wx1_1([t1:t2],1)= sparse(W)*lags([t1:t2],4);
end

for t=1:T
    t1=1+(t-1)*nc;t2=t*nc;
    Wx2_1([t1:t2],1)= sparse(W)*lags([t1:t2],8);
end

for t=1:T
    t1=1+(t-1)*nc;t2=t*nc;
    Wy1([t1:t2],1)= sparse(W)*y1([t1:t2],1);
end

% Note: There are plenty of OLS or TSLS matlab programs out there (see Lesage's
% collection of code, for example). I include my own code *in* the program
% for my own satisfaction (as a way to keep track of things).  
% You could easily cut lines 184 through 270 and go straigt to the loop
% below (which uses Lesage's hwhite(y,X) program).  

% 3) OLS and TSLS
% First State estimation
% [n,m] = size(Wy);
% cc = ones(n,1);
 %WX = [cc Wx1 Wx2 Wx4 Wy1 lags y1 x1 x2 x3 x5 x4 x6 x7];
 %[n,k] = size(WX);
 %bhat = inv(WX'*WX)*(WX'*Wy);  
 %Wybar = mean(Wy);
 %e_fs = Wy - WX*bhat;  
 %Wy_Instr = WX*bhat;
 %sse = e_fs'*e_fs/(n-k);  
 %varco = sse*inv(WX'*WX);  
 %se = sqrt(diag(varco));
 %tstat = bhat./se;
 %pvalue = 2*(1 - tcdf(abs(tstat),n-k));  
% FirstStage_results.Wy_Instr = Wy_Instr;
% FirstStage_results.bhat = bhat;
 %FirstStage_results.e_fs = e_fs;
% FirstStage_results.se = se;
 %FirstStage_results.tstat = tstat;
% FirstStage_results.pvalue = pvalue;
% WyI = Wy_Instr;

%Second Stage --------------------------------------------
% [n,m] = size(y);
% cc = ones(n,1);
% X = [cc y1 Wy x1 x2 x3 x5 x4 x6 x7];
% [n,k] = size(X);
% bhat = inv(X'*X)*(X'*y);  
% ybar = mean(y);
% e = y - X*bhat;  
% y_pred = X*bhat;
% sse = e'*e/(n-k);  
% varco = sse*inv(X'*X);  
% se = sqrt(diag(varco));
% tstat = bhat./se;
% pvalue = 2*(1 - tcdf(abs(tstat),n-k));  

%__________________________________________________________________________
%White's Corrected Standard Errors (from Lesage's Toolbox)
% white_se = hwhite(y,X);
% white_tstat = white_se.tstat;

% Much of this code may be redundant or useless to some; I include it to
% provide a nice "printout" in the command window.

% results.k = k;    
% results.n = n;
% results.X = X;
% results.y_pred = y_pred;
% results.sse = sse;
% results.bhat = bhat;
% results.rho = bhat(3,:);
% results.e = e;
% results.se = se;
% results.tstat = tstat;
% results.pvalue = pvalue;
%calclute the white se's
% se_white = bhat./white_tstat;
% pvalue_white = 2*(1 - tcdf(abs(white_tstat),n-k));
%end calculation for those stats
%  fprintf('   ***********Coefficient estimates****************\n');
%  fprintf('   ***************************************************\n');
%  fprintf('   BETA    Robust S.Errors    Robust T-stat    p-value\n');
%  B=results.bhat;
%  for i=1:length(results.bhat)
%      fprintf(' %10.3f    %10.3f     %10.3f    %10.4f\n',B(i),se_white(i),white_tstat(i),pvalue_white(i));
%  end
%  fprintf('   ***************************************************\n');
%_______________Regression Statistics__________________
% rsq = ((bhat'*X'*y)-(n*ybar^2))/(y'*y-(n*ybar^2)); %R2: Goodness of Fit
% r2adj = 1-(1-rsq)*((n-1)/(n-k));
% F = (rsq/(k-1))/((1-rsq)/(n-k));  %F-test, test of statistical significance of R2
% pvalue_F = 1-fcdf(abs(F),k-1,n-k);
% results.rsq = rsq;
% results.r2adj = r2adj;
% results.F = F;
% results.pvalue_F = pvalue_F;
% fprintf('   Model statistics\n');
% fprintf('   No. of obs.: %.0f \n',results.n);
% fprintf('   R-squared  : %.3f \n',results.rsq);
% fprintf('   Adj. R-sqr.: %.3f \n',results.r2adj);
% fprintf('   F-Statistic: %.3f \n',results.F);
% fprintf('   pvalue_F: %.4f \n',results.pvalue_F);
% fprintf('   ***************************************************\n');
% fprintf('   ***************************************************\n');


% 4) Estimating the Linear Projections: This code is a simple loop that
%    estimates the leads of the dependent variable on the regressors.  

T=179;  %number of time periods: 90 - one lag
N=26; % number of regions
cc = ones(n,1);
nobs=N*T;
% Execute the estimation for the leads of the dependent variable
% First drop the last nlost observations from both the ylead and x matrix
% Note: The loop below is for only one regressor. Either copy the code for
% any additional regressor of interest, or just change the row number for
% the "beta" and "se" lines.  

 %x = [cc y1 WyI x1 x2 x3 x5 x4 x6 x7];%IV for spatial regressor
x = [cc y1 Wy x1 x2 x3 x5 x4 x6 x7]; %Basic OLS

% Note thw code below is written assuming that line 122 has been executed
for t=1:60 %Length of the horizon
    nt = N*(T-t);
    ylead = yleads(1:nt,t);
    xlead = x(1:nt,:);
    spat_lead(t,1) = hwhite(ylead,xlead); 
end
for t=1:60
    irf(t,1) = spat_lead(t,1).beta(3,:); %Collect the coefficient on the third term in x
    se(t,1) = spat_lead(t,1).beta(3,:)/spat_lead(t,1).tstat(3,:);
end

irf_se_plus2 = irf+(se*2);
irf_se_minus2 = irf-(se*2);
irf_spatial_withbands = [irf_se_plus2 irf irf_se_minus2];

%Send irf_spatial_withbands to an Excel file if you wish:
%xlswrite('irf_spatial.xls', irf_spatial_withbands)

%The plot produced below is simple; modify as you see fit.
% Create figure
figure1 = figure;
% Create axes
axes('Parent',figure1);
box('on');
hold('all');
% Create plot
%plot(irf_spatial_withbands);
xmarker=linspace(1,60,60);
ymarker=irf_spatial_withbands;
axis([1 60 -1 1]);
plot(xmarker,irf_se_plus2,'k',xmarker,irf_se_plus2, 'ks','LineWidth',1,'MarkerSize',5);
%plot(irf_spatial_withbands,irf_spatial_withbands,'*');
hold on;
plot(xmarker,irf,'k','LineWidth',1);
hold on;
plot(xmarker,irf_se_minus2,'k',xmarker,irf_se_minus2, 'ks','LineWidth',1,'MarkerSize',5);
hold on;
plot(xmarker,linspace(0,0,60),'LineWidth',1.5);






