%
% unemployment rate regression
% including standard errors of results reported in paper
% including results of a series of test statistics
%
A=xlsread('x:\Leerstoel Regionale Economie\Unemployment.xlsx');
W1=xlsread('x:\Leerstoel Regionale Economie\BCNL.xlsx');
W=normw(W1); %normw is a self-made matlab program which should be posted in the default directory of matlab
N=12; % number of spatial units
T=41; % number of time periods
nobs=N*T;
fid=1;
y=A(:,2);%dependent variable
ym=zeros(T,1);
x=zeros(nobs,N);
for t=1:T
    t1=(t-1)*N+1;t2=t*N;
    for i=1:N
        xc(t1-1+i,i)=1;
        x(t1-1+i,i)=A(t1-1+i,3);
    end
    ym(t)=mean(y(t1:t2));    
end
for t=1:T
    t1=(t-1)*N+1;t2=t*N;
    Wy(t1:t2,1)=W*y(t1:t2,1);
end
%
% Simultaneous approach
% Produces results in Table 4 of RSUE 2016 paper
%
info.model=1; % spatial fixed efects
info.fe=0;% do not print fixed effects and their t-values
resultssar = sar_panel_FE(y(N+1:end),[y(1:end-N) Wy(1:end-N) x(N+1:end,:) x(1:end-N,:)],W,T-1,info);
vnames=strvcat('unemployment','ureg[-1]','W*ureg[-1]',...
   'unat1','unat2','unat3','unat4','unat5','unat6','unat7','unat8','unat9','unat10','unat11','unat12',...
   'u-1','u-2','u-3','u-4','u-5','u-6','u-7','u-8','u-9','u-10','u-11','u-12'); % specify variable names
%prt_sp(resultssar,vnames,1)
logliksimsfe=resultssar.lik;
info.tl=1;
info.stl=1;
info.ted=1; 
info.dyn=1;
results1=sar_jihai(y(1:end),[x(N+1:end,:) x(1:end-N,:)],W,info);
resultssar.beta=results1.theta1(1:end-2);
resultssar.rho=results1.theta1(end-1);
resultssar.tstat=results1.tstat(1:end-1);
prt_sp(resultssar,vnames,1);
[results1.theta1 results1.theta1./results1.tstat]
%
%
% Two-step approach
% Produces results in column (3) of Table 2 of RSUE 2016 paper
%
% First step
xconstant=ones(N*T,1);
[nobs K]=size([xc x]);
results=ols(y,[xc x]);
yres=results.resid;
vnames=strvcat('regional unemployment','con1','con2','con3','con4','con5','con6','con7','con8','con9','con10','con11','con12',...
'unat1','unat2','unat3','unat4','unat5','unat6','unat7','unat8','unat9','unat10','unat11','unat12'); % specify variable names
prt_reg(results,vnames);
sige=results.sige*((nobs-K)/nobs);
loglikols=-nobs/2*log(2*pi*sige)-1/(2*sige)*results.resid'*results.resid
yres=results.resid;
gam0r=results.beta(1:N);
gamr=results.beta(N+1:2*N);
%
% Second step
%
for t=1:T
    t1=(t-1)*N+1;t2=t*N;
    Wyres(t1:t2,1)=W*yres(t1:t2,1);
    Wy(t1:t2,1)=W*y(t1:t2,1);
end
info.lflag=0;
info.model=0;
resultssar = sar_panel_FE(yres(N+1:end),[ones(N*(T-1),1) yres(1:end-N) Wyres(1:end-N)],W,T-1,info);
vnames=strvcat('unemployment','constant','ureg[-1]','W*ureg[-1]'); % specify variable names
prt_sp(resultssar,vnames,1)