六步学会用MATLAB做空间计量回归详细步骤.doc

实用标准文档
1.excel与MATLAB链接：
Excel：
选项——加载项——COM加载项——转到——没有勾选项

2. MATLAB安装目录中寻找toolbox——exlink——点击,启用宏
E:\MATLAB\toolbox\exlink

然后，Excel中就出现MATLAB工具

（注意Excel中的数据：）

3.启动matlab

（1） 点击start MATLAB
（2） senddata to matlab ，并对变量矩阵变量进行命名（注意：选取变量为数值，不包括各变量）

(data表中数据进行命名)

（空间权重进行命名）
（3） 导入MATLAB中的两个矩阵变量就可以看见

4.将elhorst和jplv7两个程序文件夹复制到MATLAB安装目录的toolbox文件夹

5.设置路径：

6.输入程序，得出结果

T=30;
N=46;
W=normw(W1);
y=A(:,3);
x=A(:,[4,6]);
xconstant=ones(N*T,1);
[nobs K]=size(x);
results=ols(y,[xconstant x]);
vnames=strvcat('logcit','intercept','logp','logy');
prt_reg(results,vnames,1);
sige=results.sige*((nobs-K)/nobs);
loglikols=-nobs/2*log(2*pi*sige)-1/(2*sige)*results.resid'*results.resid
% The (robust)LM tests developed by Elhorst
LMsarsem_panel(results,W,y,[xconstant x]); % (Robust) LM tests
解释
每一行分别表示：
该面板数据的时期数为30（T=30），
该面板数据有30个地区（N=30），
将空间权重矩阵标准化（W=normw(w1)），
将名为A（以矩阵形式出现在MATLABA中）的变量的第3列数据定义为被解释变量y,
将名为A的变量的第4、5、6列数据定义为解释变量矩阵x，
定义一个有N*T行，1列的全1矩阵，该矩阵名为：xconstant，（ones即为全1矩阵）
说明解释变量矩阵x的大小：有nobs行，K列。（size为描述矩阵的大小）。

附录：
静态面板空间计量经济学

一、OLS静态面板编程
1、普通面板编程
T=30;
N=46;
W=normw(W1);
y=A(:,3);
x=A(:,[4,6]);
xconstant=ones(N*T,1);
[nobs K]=size(x);
results=ols(y,[xconstant x]);
vnames=strvcat('logcit','intercept','logp','logy');
prt_reg(results,vnames,1);
sige=results.sige*((nobs-K)/nobs);
loglikols=-nobs/2*log(2*pi*sige)-1/(2*sige)*results.resid'*results.resid
% The (robust)LM tests developed by Elhorst
LMsarsem_panel(results,W,y,[xconstant x]); % (Robust) LM tests

2、空间固定OLS （spatial-fixed effects）
T=30;
N=46;
W=normw(W1);
y=A(:,3);
x=A(:,[4,6]);
xconstant=ones(N*T,1);
[nobs K]=size(x);
model=1;
[ywith,xwith,meanny,meannx,meanty,meantx]=demean(y,x,N,T,model);
results=ols(ywith,xwith);
vnames=strvcat('logcit','logp','logy'); % should be changed if x is changed
prt_reg(results,vnames);
sfe=meanny-meannx*results.beta; % including the constant term
yme = y - mean(y);
et=ones(T,1);
error=y-kron(et,sfe)-x*results.beta;
rsqr1 = error'*error;
rsqr2 = yme'*yme;
FE_rsqr2 = 1.0 - rsqr1/rsqr2 % r-squared including fixed effects
sige=results.sige*((nobs-K)/nobs);
logliksfe=-nobs/2*log(2*pi*sige)-1/(2*sige)*results.resid'*results.resid
LMsarsem_panel(results,W,ywith,xwith); % (Robust) LM tests

3、时期固定OLS（time-period fixed effects）
T=30;
N=46;
W=normw(W1);
y=A(:,3);
x=A(:,[4,6]);
xconstant=ones(N*T,1);
[nobs K]=size(x);
model=2;
[ywith,xwith,meanny,meannx,meanty,meantx]=demean(y,x,N,T,model);
results=ols(ywith,xwith);
vnames=strvcat('logcit','logp','logy'); % should be changed if x is changed
prt_reg(results,vnames);
tfe=meanty-meantx*results.beta; % including the constant term
yme = y - mean(y);
en=ones(N,1);
error=y-kron(tfe,en)-x*results.beta;
rsqr1 = error'*error;
rsqr2 = yme'*yme;
FE_rsqr2 = 1.0 - rsqr1/rsqr2 % r-squared including fixed effects
sige=results.sige*((nobs-K)/nobs);
logliktfe=-nobs/2*log(2*pi*sige)-1/(2*sige)*results.resid'*results.resid
LMsarsem_panel(results,W,ywith,xwith); % (Robust) LM tests

4、空间与时间双固定模型
T=30;
N=46;
W=normw(W1);
y=A(:,3);
x=A(:,[4,6]);
xconstant=ones(N*T,1);
[nobs K]=size(x);
model=3;
[ywith,xwith,meanny,meannx,meanty,meantx]=demean(y,x,N,T,model);
results=ols(ywith,xwith);
vnames=strvcat('logcit','logp','logy'); % should be changed if x is changed
prt_reg(results,vnames)
en=ones(N,1);
et=ones(T,1);
intercept=mean(y)-mean(x)*results.beta;
sfe=meanny-meannx*results.beta-kron(en,intercept);
tfe=meanty-meantx*results.beta-kron(et,intercept);
yme = y - mean(y);
ent=ones(N*T,1);
error=y-kron(tfe,en)-kron(et,sfe)-x*results.beta-kron(ent,intercept);
rsqr1 = error'*error;
rsqr2 = yme'*yme;
FE_rsqr2 = 1.0 - rsqr1/rsqr2 % r-squared including fixed effects
sige=results.sige*((nobs-K)/nobs);
loglikstfe=-nobs/2*log(2*pi*sige)-1/(2*sige)*results.resid'*results.resid
LMsarsem_panel(results,W,ywith,xwith); % (Robust) LM tests

二、静态面板SAR模型
1、无固定效应（No fixed effects）
T=30;
N=46;
W=normw(W1);
y=A(:,[3]);
x=A(:,[4,6]);
for t=1:T
t1=(t-1)*N+1;t2=t*N;
wx(t1:t2,:)=W*x(t1:t2,:);
end
xconstant=ones(N*T,1);
[nobs K]=size(x);
info.lflag=0;
info.model=0;
info.fe=0;
results=sar_panel_FE(y,[xconstant x],W,T,info);
vnames=strvcat('logcit','intercept','logp','logy');
prt_spnew(results,vnames,1)

% Print out effects estimates
spat_model=0;
direct_indirect_effects_estimates(results,W,spat_model);
panel_effects_sar(results,vnames,W);
2、空间固定效应（Spatial fixed effects）
T=30;
N=46;
W=normw(W1);
y=A(:,[3]);
x=A(:,[4,6]);
for t=1:T
t1=(t-1)*N+1;t2=t*N;
wx(t1:t2,:)=W*x(t1:t2,:);
end
xconstant=ones(N*T,1);
[nobs K]=size(x);
info.lflag=0;
info.model=1;
info.fe=0;
results=sar_panel_FE(y,x,W,T,info);
vnames=strvcat('logcit','logp','logy');
prt_spnew(results,vnames,1)

% Print out effects estimates
spat_model=0;
direct_indirect_effects_estimates(results,W,spat_model);
panel_effects_sar(results,vnames,W);
3、时点固定效应（Time period fixed effects）
T=30;
N=46;
W=normw(W1);
y=A(:,[3]);
x=A(:,[4,6]);
for t=1:T
t1=(t-1)*N+1;t2=t*N;
wx(t1:t2,:)=W*x(t1:t2,:);
end
xconstant=ones(N*T,1);
[nobs K]=size(x);
info.lflag=0; % required for exact results
info.model=2;
info.fe=0; % Do not print intercept and fixed effects; use info.fe=1 to turn on
results=sar_panel_FE(y,x,W,T,info);
vnames=strvcat('logcit','logp','logy');
prt_spnew(results,vnames,1)

% Print out effects estimates
spat_model=0;
direct_indirect_effects_estimates(results,W,spat_model);
panel_effects_sar(results,vnames,W);
4、双固定效应（Spatial and time period fixed effects）
T=30;
N=46;
W=normw(W1);
y=A(:,[3]);
x=A(:,[4,6]);
for t=1:T
t1=(t-1)*N+1;t2=t*N;
wx(t1:t2,:)=W*x(t1:t2,:);
end
xconstant=ones(N*T,1);
[nobs K]=size(x);
info.lflag=0; % required for exact results
info.model=3;
info.fe=0; % Do not print intercept and fixed effects; use info.fe=1 to turn on
results=sar_panel_FE(y,x,W,T,info);
vnames=strvcat('logcit','logp','logy');
prt_spnew(results,vnames,1)
% Print out effects estimates
spat_model=0;
direct_indirect_effects_estimates(results,W,spat_model);
panel_effects_sar(results,vnames,W);

三、静态面板SDM模型
1、无固定效应（No fixed effects）
T=30;
N=46;
W=normw(W1);
y=A(:,[3]);
x=A(:,[4,6]);
for t=1:T
t1=(t-1)*N+1;t2=t*N;
wx(t1:t2,:)=W*x(t1:t2,:);
end
xconstant=ones(N*T,1);
[nobs K]=size(x);
info.lflag=0;
info.model=0;
info.fe=0;
results=sar_panel_FE(y,[xconstant x wx],W,T,info);
vnames=strvcat('logcit','intercept','logp','logy','W*logp','W*logy');
prt_spnew(results,vnames,1)
% Print out effects estimates
spat_model=1;
direct_indirect_effects_estimates(results,W,spat_model);
panel_effects_sdm(results,vnames,W);
2、空间固定效应（Spatial fixed effects）
T=30;
N=46;
W=normw(W1);
y=A(:,[3]);
x=A(:,[4,6]);
for t=1:T
t1=(t-1)*N+1;t2=t*N;
wx(t1:t2,:)=W*x(t1:t2,:);
end
xconstant=ones(N*T,1);
[nobs K]=size(x);
info.lflag=0; % required for exact results
info.model=1;
info.fe=0; % Do not print intercept and fixed effects; use info.fe=1 to turn on
results=sar_panel_FE(y,[x wx],W,T,info);
vnames=strvcat('logcit','logp','logy','W*logp','W*logy');
prt_spnew(results,vnames,1)
% Print out effects estimates
spat_model=1;
direct_indirect_effects_estimates(results,W,spat_model);
panel_effects_sdm(results,vnames,W);
3、时点固定效应（Time period fixed effects）
T=30;
N=46;
W=normw(W1);
y=A(:,[3]);
x=A(:,[4,6]);
for t=1:T
t1=(t-1)*N+1;t2=t*N;
wx(t1:t2,:)=W*x(t1:t2,:);
end
xconstant=ones(N*T,1);
[nobs K]=size(x);
info.lflag=0; % required for exact results
info.model=2;
info.fe=0; % Do not print intercept and fixed effects; use info.fe=1 to turn on
% New routines to calculate effects estimates
results=sar_panel_FE(y,[x wx],W,T,info);
vnames=strvcat('logcit','logp','logy','W*logp','W*logy');
% Print out coefficient estimates
prt_spnew(results,vnames,1)
% Print out effects estimates
spat_model=1;
direct_indirect_effects_estimates(results,W,spat_model);
panel_effects_sdm(results,vnames,W)
4、双固定效应（Spatial and time period fixed effects）
T=30;
N=46;
W=normw(W1);
y=A(:,[3]);
x=A(:,[4,6]);
for t=1:T
t1=(t-1)*N+1;t2=t*N;
wx(t1:t2,:)=W*x(t1:t2,:);
end
xconstant=ones(N*T,1);
[nobs K]=size(x);
info.bc=0;
info.lflag=0; % required for exact results
info.model=3;
info.fe=0; % Do not print intercept and fixed effects; use info.fe=1 to turn on
results=sar_panel_FE(y,[x wx],W,T,info);
vnames=strvcat('logcit','logp','logy','W*logp','W*logy');
prt_spnew(results,vnames,1)
% Print out effects estimates
spat_model=1;
direct_indirect_effects_estimates(results,W,spat_model);
panel_effects_sdm(results,vnames,W)

wald test spatial lag
% Wald test for spatial Durbin model against spatial lag model
btemp=results.parm;
varcov=results.cov;
Rafg=zeros(K,2*K+2);
for k=1:K
Rafg(k,K+k)=1; % R(1,3)=0 and R(2,4)=0;
end
Wald_spatial_lag=(Rafg*btemp)'*inv(Rafg*varcov*Rafg')*Rafg*btemp
prob_spatial_lag=1-chis_cdf (Wald_spatial_lag, K)

wald test spatial error
% Wald test spatial Durbin model against spatial error model
R=zeros(K,1);
for k=1:K
R(k)=btemp(2*K+1)*btemp(k)+btemp(K+k); % k changed in 1, 7/12/2010
% R(1)=btemp(5)*btemp(1)+btemp(3);
% R(2)=btemp(5)*btemp(2)+btemp(4);
end
Rafg=zeros(K,2*K+2);
for k=1:K
Rafg(k,k) =btemp(2*K+1); % k changed in 1, 7/12/2010
Rafg(k,K+k) =1;
Rafg(k,2*K+1)=btemp(k);
% Rafg(1,1)=btemp(5);Rafg(1,3)=1;Rafg(1,5)=btemp(1);
% Rafg(2,2)=btemp(5);Rafg(2,4)=1;Rafg(2,5)=btemp(2);
end
Wald_spatial_error=R'*inv(Rafg*varcov*Rafg')*R
prob_spatial_error=1-chis_cdf (Wald_spatial_error,K)

LR test spatial lag
resultssar=sar_panel_FE(y,x,W,T,info);
LR_spatial_lag=-2*(resultssar.lik-results.lik)
prob_spatial_lag=1-chis_cdf (LR_spatial_lag,K)

LR test spatial error
resultssem=sem_panel_FE(y,x,W,T,info);
LR_spatial_error=-2*(resultssem.lik-results.lik)
prob_spatial_error=1-chis_cdf (LR_spatial_error,K)

5、空间随机效应与时点固定效应模型
T=30;
N=46;
W=normw(W1);
y=A(:,[3]);
x=A(:,[4,6]);
for t=1:T
t1=(t-1)*N+1;t2=t*N;
wx(t1:t2,:)=W*x(t1:t2,:);
end
xconstant=ones(N*T,1);
[nobs K]=size(x);
[ywith,xwith,meanny,meannx,meanty,meantx]=demean(y,[x wx],N,T,2); % 2=time dummies
info.model=1;
results=sar_panel_RE(ywith,xwith,W,T,info);
prt_spnew(results,vnames,1)

spat_model=1;
direct_indirect_effects_estimates(results,W,spat_model);
panel_effects_sdm(results,vnames,W)

wald test spatial lag
btemp=results.parm(1:2*K+2);
varcov=results.cov(1:2*K+2,1:2*K+2);
Rafg=zeros(K,2*K+2);
for k=1:K
Rafg(k,K+k)=1; % R(1,3)=0 and R(2,4)=0;
end
Wald_spatial_lag=(Rafg*btemp)'*inv(Rafg*varcov*Rafg')*Rafg*btemp
prob_spatial_lag= 1-chis_cdf (Wald_spatial_lag, K)

wald test spatial error
R=zeros(K,1);
for k=1:K
R(k)=btemp(2*K+1)*btemp(k)+btemp(K+k); % k changed in 1, 7/12/2010
% R(1)=btemp(5)*btemp(1)+btemp(3);
% R(2)=btemp(5)*btemp(2)+btemp(4);
end
Rafg=zeros(K,2*K+2);
for k=1:K
Rafg(k,k) =btemp(2*K+1); % k changed in 1, 7/12/2010
Rafg(k,K+k) =1;
Rafg(k,2*K+1)=btemp(k);
% Rafg(1,1)=btemp(5);Rafg(1,3)=1;Rafg(1,5)=btemp(1);
% Rafg(2,2)=btemp(5);Rafg(2,4)=1;Rafg(2,5)=btemp(2);
end
Wald_spatial_error=R'*inv(Rafg*varcov*Rafg')*R
prob_spatial_error= 1-chis_cdf (Wald_spatial_error,K)

LR test spatial lag
resultssar=sar_panel_RE(ywith,xwith(:,1:K),W,T,info);
LR_spatial_lag=-2*(resultssar.lik-results.lik)
prob_spatial_lag=1-chis_cdf (LR_spatial_lag,K)

LR test spatial error
resultssem=sem_panel_RE(ywith,xwith(:,1:K),W,T,info);
LR_spatial_error=-2*(resultssem.lik-results.lik)
prob_spatial_error=1-chis_cdf (LR_spatial_error,K)

四、静态面板SEM模型
1、无固定效应（No fixed effects）
T=30;
N=46;
W=normw(W1);
y=A(:,[3]);
x=A(:,[4,6]);
for t=1:T
t1=(t-1)*N+1;t2=t*N;
wx(t1:t2,:)=W*x(t1:t2,:);
end
xconstant=ones(N*T,1);
[nobs K]=size(x);
info.lflag=0;
info.model=0;
info.fe=0;
results=sem_panel_FE(y,[xconstant x],W,T,info);
vnames=strvcat('logcit','intercept','logp','logy');
prt_spnew(results,vnames,1)
% Print out effects estimates
spat_model=0;
direct_indirect_effects_estimates(results,W,spat_model);
panel_effects_sar(results,vnames,W);
2、空间固定效应（Spatial fixed effects）
T=30;
N=46;
W=normw(W1);
y=A(:,[3]);
x=A(:,[4,6]);
for t=1:T
t1=(t-1)*N+1;t2=t*N;
wx(t1:t2,:)=W*x(t1:t2,:);
end
xconstant=ones(N*T,1);
[nobs K]=size(x);
info.lflag=0;
info.model=1;
info.fe=0;
results=sem_panel_FE(y,x,W,T,info);
vnames=strvcat('logcit','logp','logy');
prt_spnew(results,vnames,1)
% Print out effects estimates
spat_model=0;
direct_indirect_effects_estimates(results,W,spat_model);
panel_effects_sar(results,vnames,W);
3、时点固定效应（Time period fixed effects）
T=30;
N=46;
W=normw(W1);
y=A(:,[3]);
x=A(:,[4,6]);
for t=1:T
t1=(t-1)*N+1;t2=t*N;
wx(t1:t2,:)=W*x(t1:t2,:);
end
xconstant=ones(N*T,1);
[nobs K]=size(x);
info.lflag=0; % required for exact results
info.model=2;
info.fe=0; % Do not print intercept and fixed effects; use info.fe=1 to turn on
results=sem_panel_FE(y,x,W,T,info);
vnames=strvcat('logcit','logp','logy');
prt_spnew(results,vnames,1)
% Print out effects estimates
spat_model=0;
direct_indirect_effects_estimates(results,W,spat_model);
panel_effects_sar(results,vnames,W);
4、双固定效应（Spatial and time period fixed effects）
T=30;
N=46;
W=normw(W1);
y=A(:,[3]);
x=A(:,[4,6]);
for t=1:T
t1=(t-1)*N+1;t2=t*N;
wx(t1:t2,:)=W*x(t1:t2,:);
end
xconstant=ones(N*T,1);
[nobs K]=size(x);
info.lflag=0; % required for exact results
info.model=3;
info.fe=0; % Do not print intercept and fixed effects; use info.fe=1 to turn on
results=sem_panel_FE(y,x,W,T,info);
vnames=strvcat('logcit','logp','logy');
prt_spnew(results,vnames,1)
% Print out effects estimates
spat_model=0;
direct_indirect_effects_estimates(results,W,spat_model);
panel_effects_sar(results,vnames,W);

五、静态面板SDEM模型
1、无固定效应（No fixed effects）
T=30;
N=46;
W=normw(W1);
y=A(:,[3]);
x=A(:,[4,6]);
for t=1:T
t1=(t-1)*N+1;t2=t*N;
wx(t1:t2,:)=W*x(t1:t2,:);
end
xconstant=ones(N*T,1);
[nobs K]=size(x);
info.lflag=0;
info.model=0;
info.fe=0;
results=sem_panel_FE(y,[xconstant x wx],W,T,info);
vnames=strvcat('logcit','intercept','logp','logy','W*logp','W*logy');
prt_spnew(results,vnames,1)
% Print out effects estimates
spat_model=1;
direct_indirect_effects_estimates(results,W,spat_model);
panel_effects_sdm(results,vnames,W);
2、空间固定效应（Spatial fixed effects）
T=30;
N=46;
W=normw(W1);
y=A(:,[3]);
x=A(:,[4,6]);
for t=1:T
t1=(t-1)*N+1;t2=t*N;
wx(t1:t2,:)=W*x(t1:t2,:);
end
xconstant=ones(N*T,1);
[nobs K]=size(x);
info.lflag=0; % required for exact results
info.model=1;
info.fe=0; % Do not print intercept and fixed effects; use info.fe=1 to turn on
results=sem_panel_FE(y,[x wx],W,T,info);
vnames=strv