%South CA n=8 in Brady (2011)

load South
% load Sample1_0 %diff. from Sample1 in that 0.0012=current housing price

%my prepared data loaded from the Matlab workspace (previously saved as a .mat file)

T = 56;  %Time in months for each region
nc = 31; %The number of regions

reg = data(:,26); %County numbers
dt = data(:,2); %date
y = data(:,3); %housing prices (natural log)
y1 = data(:,4); %housing prices lag
x1 = data(:,8);  %unemployment
x2 = data(:,11);  %construction (natural log)
x3 = data(:,14);  %realrate
x4 = data(:,18);  %population (natural log)
x5 = data(:,22); %industrial production (natural log)
x6 = data(:,27:29);  %season dummies
x7 = data(:,38:67);  %county dummies (22 for 23 counties)

%lags of the regressors
y2 = data(:,5); y3 = data(:,6); y4 = data(:,7); %more housing lags
x1_1 = data(:,9); x1_2 = data(:,10);%unemployment lags
x2_1 = data(:,12); x2_2 = data(:,13); %construction lags
x3_1 = data(:,15); x3_2 = data(:,16);%realrate lags
x4_1 = data(:,19); x4_2 = data(:,20);%population lags
x5_1 = data(:,23); x5_2 = data(:,24); %IP lags

%Housing price leads +1 through +36 (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."

W=[1,1,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,1,1,0,1,0,1,0,0,0,1,0,0,0;1,1,0,0,0,0,0,1,0,0,1,0,0,0,0,1,0,0,1,1,0,0,0,0,0,1,1,0,0,0,0;0,0,1,0,0,0,0,0,0,0,0,0,0,1,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0;0,0,0,1,0,0,1,0,1,1,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,1,0,0;0,0,0,0,1,1,0,0,0,1,0,0,0,0,0,0,0,0,0,0,1,0,1,0,0,0,0,0,1,1,0;0,0,0,0,1,1,0,0,0,0,0,0,1,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0;0,0,0,1,0,0,1,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0;1,1,0,0,0,0,0,1,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,1,1,0,0,0,0;0,0,0,1,0,0,1,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,1,0,0,0,1,0,0,0;0,0,0,1,1,0,0,0,0,1,0,0,0,0,0,0,1,0,0,0,1,0,0,1,1,0,0,0,0,0,0;0,1,0,0,0,0,0,1,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,0,0,0,1;0,0,0,0,0,0,0,0,0,0,0,1,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0;0,0,0,0,0,1,0,0,0,0,0,0,1,0,1,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0;0,0,1,0,0,0,0,0,0,0,0,1,0,1,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0;0,0,0,0,0,1,0,0,0,0,0,0,1,0,1,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0;0,1,1,0,0,0,0,0,0,0,0,0,0,1,0,1,0,0,0,1,0,0,0,0,0,1,0,0,0,0,1;0,0,0,1,0,0,0,0,1,1,0,0,0,0,0,0,1,0,0,0,0,0,0,1,1,0,0,0,0,0,0;0,0,0,0,0,0,0,0,0,0,0,0,1,0,1,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0;1,1,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,1,0,0,1,0,0,0,0,0,0,0,0,0;1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,1,0,0,0,1,0,0,0,1,0,0,0;0,0,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,1,0,1,0,0,0,0,0,0,0,0;1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,1,0,1,1,0,0,0,0,0,0;0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,1,0,0,0,0,0,0,1,0;1,0,0,0,0,0,0,0,1,1,0,0,0,0,0,0,1,0,0,1,0,1,0,1,1,0,0,1,0,0,0;0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,1,0,0,0,0,1,0,1,1,0,0,0,0,0,0;0,1,0,0,0,0,0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,0,0,1,1,0,0,0,1;0,1,0,0,0,0,0,1,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,0,0,0,0;1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,1,0,0,0,1,0,0,0,1,0,0,0;0,0,0,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0;0,0,0,0,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,1,0;0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,0,0,1,0,0,0,0,1;];

x = zeros(31,31);

for i=1:31
    for j=1:31
        if j~=i
    x(i,j)=W(i,j);
        end
    end
end


% In the next few lines I sort the data
lags = [y2 y3 y4 x1_1  x1_2 x2_1 x2_2 x3_1 x3_2 x4_1 x4_2 x5_1 x5_2];
all = [dt reg y y1 x1 x2 x3 x4 x5 x6 x7 lags];
[n,m] = size(all);

% 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) . . . )


% I get rid of 199804 obs. since for county 21, the laged housing price is
% imputed as 0.0012.

all = sortrows(all,1); %This sorts the matrix according to data ("dt"), this is required by Lesage's code used below.
all(1:31,:)=[]; % delete the obs in 199804.

[n,m] = size(all); 
y =  (all(:,3))*100; %Here I just redefine the variables in the "all" matrix
y1 = (all(:,4))*100;
x1 = all(:,5);
x2 = (all(:,6))*100;
x3 = all(:,7);
x4 = (all(:,8))*100;
x5 = (all(:,9))*100;
x6 = all(:,10:12);
x7 = all(:,13:42);
dt = all(:,1);
lags = (all(:,43:55))*100;
W = x; %Rename the weighting matrix consistent with Lesage's code.

y2=lags(:,1);
y3=lags(:,2);
y4=lags(:,3);
x1_1=lags(:,4);
x1_2=lags(:,5);
x2_1=lags(:,6);
x2_2=lags(:,7);
x3_1=lags(:,8);
x3_2=lags(:,9);
x4_1=lags(:,10);
x4_2=lags(:,11);
x5_1=lags(:,12);
x5_2=lags(:,13);

% Lesage's code: "sparse" "normw" and "slag" are from Lesage's spatial toolbox (see his webiste).  
WW = sparse(W);
W1 = normw(WW);
W = slag(W1,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.

% I delete the obs in 199804.

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: 

% define nt-by-nt weithting matrix

W_nt=kron(eye(89),W);


% OLS with spatial regressor Shulin Shen write the OLS code
[n,m] = size(y);
cc = ones(n,1);
X = [cc Wy y1 x1 x2 x3 x4 x5 x6 x7];
[n,k] = size(X);
bhat = inv(X'*X)*(X'*y);  
ybar = mean(y);
sst=(y-ybar)'*(y-ybar); % Shulin Shen add the total sum of squares
e_fs = y - X*bhat;  
y_OLS = X*bhat;
sse = e_fs'*e_fs/(n-k);  
varco = sse*inv(X'*X);  
se = sqrt(diag(varco));
tstat = bhat./se;
pvalue = 2*(1 - tcdf(abs(tstat),n-k));  
OLS_results.y_OLS = y_OLS;
OLS_results.bhat = bhat;
OLS_results.e_fs = e_fs;
OLS_results.se = se;
OLS_results.tstat = tstat;
OLS_results.pvalue = pvalue;
OLS_results.r2=1-(n-k)*sse/sst; % Shulin Shen add the R-square
OLS_results.ad_r2=1-(n-1)/(n-k-2)*(n-k)*sse/sst; % Shulin Shen add the adjusted R-square

% Shulin Shen's White's Corrected Standard Errors

sigma=diag(e_fs.^2);

varco_white_beta=inv(X'*X)*(X'*sigma*X)*inv(X'*X);

se_white=sqrt(diag(varco_white_beta));

tstat_white=bhat./se_white;

OLS_results.pvalue_white=2*(1 - tcdf(abs(tstat_white),n-k));

OLS_RE=[OLS_results.bhat se_white OLS_results.pvalue_white OLS_results.pvalue OLS_results.tstat OLS_results.se];


%White's Corrected Standard Errors (from Lesage's Toolbox)
% white_se = hwhite(y,X);
% white_tstat = white_se.tstat;
%calclute the white se's
% se_white = bhat./white_tstat;
% pvalue_white = 2*(1 - tcdf(abs(white_tstat),n-k));
% OLS_results.se_white = se_white;
% OLS_results.pvalue_white= pvalue_white;
% end of OLS 


% 4) IV 

W_nt=kron(eye(89),W);
W1_nt=kron(eye(89),W1);

% First Stage estimation
[n,m] = size(Wy);
cc = ones(n,1);
WX = [cc Wx1 Wx2 Wx4 Wy1 y1 x1 x2 x3 x5 x4 x6 x7]; %Shulin Shen excluded lags from the author's code
[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)); 

% First stage F-test

RSSR=(Wy-Wybar)'*(Wy-Wybar); % Restricted sum of squared residuals
USSR=e_fs'*e_fs;
F=((RSSR-USSR)/(k-1))/(USSR/(n-k)); 
F_pvalue=1-fcdf(F,k-1,n-k);
FirstStage_results.F_test=F;
FirstStage_results.F_test_pvalue=F_pvalue;
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;

First_IV=[bhat pvalue tstat se];


%Second Stage --------------------------------------------
[n,m] = size(y);
cc = ones(n,1);
X = [cc WyI y1 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;
resluts.pvalue_white = 2*(1 - tcdf(abs(white_tstat),n-k));


IV_RE=[results.bhat se_white resluts.pvalue_white results.pvalue results.tstat results.se];

% 5) IV 2

W_nt=kron(eye(89),W);
W1_nt=kron(eye(89),W1);

% First Stage estimation
[n,m] = size(Wy);
cc = ones(n,1);

x1_1=lags(:,4);
x2_1=lags(:,8);

WX = [cc Wx1 Wx2 Wx4 Wy1 Wx1_1 Wx2_1 y1 x3 x5 x4 x6 x7]; 
% WX = [cc Wx1 Wx2 Wx4 Wy1 lags y1 x3 x5 x4 x6 x7]; %Shulin Shen excluded lags from the author's code
[n,k] = size(WX);
bhat_wy = inv(WX'*WX)*(WX'*Wy); 
Wy_Instr = WX*bhat_wy;
Wybar = mean(Wy);
e_fs = Wy - WX*bhat_wy;  
sse = e_fs'*e_fs/(n-k);  
varco = sse*inv(WX'*WX);  
se_wy= sqrt(diag(varco));
tstat_wy= bhat_wy./se_wy;
pvalue_wy= 2*(1 - tcdf(abs(tstat_wy),n-k)); 
RSSR=(Wy-Wybar)'*(Wy-Wybar); % Restricted sum of squared residuals
USSR=e_fs'*e_fs;
F_wy=((RSSR-USSR)/(k-1))/(USSR/(n-k)); 
F_wy_pvalue=1-fcdf(F_wy,k-1,n-k);

First_IV_Wy=[bhat_wy pvalue_wy tstat_wy se_wy];

[n,m] = size(x1);
cc = ones(n,1);
WX = [cc Wx1 Wx2 Wx4 Wy1 Wx1_1 Wx2_1 y1 x3 x5 x4 x6 x7]; 
% WX = [cc Wx1 Wx2 Wx4 Wy1 x1_1 x2_1 y1 x3 x5 x4 x6 x7]; %Shulin Shen excluded lags from the author's code
[n,k] = size(WX);
bhat_x1=inv(WX'*WX)*(WX'*x1);
x1_Instr = WX*bhat_x1;
x1bar = mean(x1);
e_fs = x1 - WX*bhat_x1;  
sse = e_fs'*e_fs/(n-k);  
varco = sse*inv(WX'*WX);  
se_x1= sqrt(diag(varco));
tstat_x1= bhat_x1./se_x1;
pvalue_x1= 2*(1 - tcdf(abs(tstat_x1),n-k)); 
RSSR=(x1-x1bar)'*(x1-x1bar); % Restricted sum of squared residuals
USSR=e_fs'*e_fs;
F_x1=((RSSR-USSR)/(k-1))/(USSR/(n-k)); 
F_x1_pvalue=1-fcdf(F_x1,k-1,n-k);

First_IV_x1=[bhat_x1 pvalue_x1 tstat_x1 se_x1];

[n,m] = size(x2);
cc = ones(n,1);
WX = [cc Wx1 Wx2 Wx4 Wy1 Wx1_1 Wx2_1 y1 x3 x5 x4 x6 x7];
% WX = [cc Wx1 Wx2 Wx4 Wy1 x1_1 x2_1 y1 x3 x5 x4 x6 x7]; %Shulin Shen excluded lags from the author's code
[n,k] = size(WX);
bhat_x2=inv(WX'*WX)*(WX'*x2); 
x2_Instr = WX*bhat_x2;
x2bar = mean(x2);
e_fs = x2 - WX*bhat_x2;  
sse = e_fs'*e_fs/(n-k);  
varco = sse*inv(WX'*WX);  
se_x2= sqrt(diag(varco));
tstat_x2= bhat_x2./se_x2;
pvalue_x2= 2*(1 - tcdf(abs(tstat_x2),n-k)); 
RSSR=(x2-x2bar)'*(x2-x2bar); % Restricted sum of squared residuals
USSR=e_fs'*e_fs;
F_x2=((RSSR-USSR)/(k-1))/(USSR/(n-k)); 
F_x2_pvalue=1-fcdf(F_x2,k-1,n-k);
 
First_IV_x2=[bhat_x2 pvalue_x2 tstat_x2 se_x2];

%Second Stage --------------------------------------------
[n,m] = size(y);
cc = ones(n,1);
X = [cc Wy_Instr y1 x1_Instr x2_Instr 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;
resluts.pvalue_white = 2*(1 - tcdf(abs(white_tstat),n-k));


IV2_RE=[results.bhat se_white resluts.pvalue_white results.pvalue results.tstat results.se];





