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

1、1.excel与 MATLAB链接：Excel ：选项 加载项 COM 加载项 转到 没有勾 选项2. MATLAB安装目录中寻找 toolbox exlink 点击 , 启用宏E:MATLABtoolboxexlink然后，Excel 中就出 现 MATLAB工具1（注意Excel 中的数据：）3. 启动 matlab（1（点击 start MATLAB（2（senddata to matlab ，并对变 量矩 阵变 量进行命名（注意：选取变量为数值，不包括各 变量）2(data 表中数据 进行命名 )（空间权 重进行命名）（3（导入 MATLAB中的两个矩 阵变 量就可以看 见34.将 e

2、lhorst 和 jplv7 两个程序文件 夹复制到 MATLAB 安装目录的 toolbox 文件夹5.设置路径：46.输入程序，得出结果T=30;N=46;W=normw(W1);y=A(:,3);5x=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*lo

3、g(2*pi*sige)-1/(2*sige)*results.resid*results.resid% The (robust)LM tests developed by ElhorstLMsarsem_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

4、、6 列数据定 义为解释变 量矩 阵 x，定义一个有 N*T 行，1 列的全 1 矩阵，该矩阵名为：xconstant，（ones 即为全 1 矩阵）6说明解 释变量矩 阵 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

5、,logy);prt_reg(results,vnames,1);7sige=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 ElhorstLMsarsem_panel(results,W,y,xconstant x); % (Robust) LM tests2、空间固定 OLS （spatial-fixed effects）T=30;N=46;W=normw(W1);y=A

6、(:,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);8vnames=strvcat(logcit,logp,logy); % should be changed if x is changedprt_reg(results,vnames);sfe=meanny-meannx*results.beta; % including the constant

7、termyme = 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 effectssige=results.sige*(nobs-K)/nobs);logliksfe=-nobs/2*log(2*pi*sige)-1/(2*sige)*results.resid*results.residLMsarsem_panel(results,W,

8、ywith,xwith); % (Robust) LM tests3、时期固定 OLS（time-period fixed effects）T=30;N=46;9W=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 i

9、f x is changedprt_reg(results,vnames);tfe=meanty-meantx*results.beta; % including the constant termyme = y - mean(y);en=ones(N,1);error=y-kron(tfe,en)-x*results.beta;rsqr1 = error*error;rsqr2 = yme*yme;10FE_rsqr2 = 1.0 - rsqr1/rsqr2 % r-squared including fixed effectssige=results.sige*(nobs-K)/nobs)

10、;logliktfe=-nobs/2*log(2*pi*sige)-1/(2*sige)*results.resid*results.residLMsarsem_panel(results,W,ywith,xwith); % (Robust) LM tests4、空间与时间双固定模型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=o

11、ls(ywith,xwith);11vnames=strvcat(logcit,logp,logy); % should be changed if x is changedprt_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=one

12、s(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 effectssige=results.sige*(nobs-K)/nobs);loglikstfe=-nobs/2*log(2*pi*sige)-1/(2*sige)*results.resid*results.residLMsarsem_panel(re

13、sults,W,ywith,xwith); % (Robust) LM tests12二、静态面板 SAR 模型1、无固定效应（No fixed effects）T=30;N=46;W=normw(W1);y=A(:,3);x=A(:,4,6);fort=1:Tt1=(t-1)*N+1;t2=t*N;wx(t1:t2,:)=W*x(t1:t2,:);endxconstant=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);vna

14、mes=strvcat(logcit, intercept, logp, logy);prt_spnew(results,vnames,1)% Print out effects estimates spat_model=0;13direct_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);fort=1:Tt1=(t-1)*N+1

15、;t2=t*N;wx(t1:t2,:)=W*x(t1:t2,:);endxconstant=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_mo

16、del);14panel_effects_sar(results,vnames,W);3、时点固定效 应（Time period fixed effects）T=30;N=46;W=normw(W1);y=A(:,3);x=A(:,4,6);fort=1:Tt1=(t-1)*N+1;t2=t*N;wx(t1:t2,:)=W*x(t1:t2,:);endxconstant=ones(N*T,1);nobs K=size(x);info.lflag=0;% required for exact resultsinfo.model=2;info.fe=0;% Do not print interce

17、pt and fixed effects; use info.fe=1to turn onresults=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);15panel_effects_sar(results,vnames,W);4、双固定效应（Spatial and t

18、ime period fixed effects）T=30;N=46;W=normw(W1);y=A(:,3);x=A(:,4,6);fort=1:Tt1=(t-1)*N+1;t2=t*N;wx(t1:t2,:)=W*x(t1:t2,:);endxconstant=ones(N*T,1);nobs K=size(x);info.lflag=0;% required for exact resultsinfo.model=3;info.fe=0;% Do not print intercept and fixed effects; use info.fe=1to turn onresults=s

19、ar_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);16三、静态面板 SDM 模型1、无固定效应（No fixed effects）T=30;N=46;W=normw(W1);y=A(:,3);x=A(:

20、,4,6);fort=1:Tt1=(t-1)*N+1;t2=t*N;wx(t1:t2,:)=W*x(t1:t2,:);endxconstant=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 estimate

21、s spat_model=1;direct_indirect_effects_estimates(results,W,spat_model);17panel_effects_sdm(results,vnames,W);2、空间固定效应（Spatial fixed effects）T=30;N=46;W=normw(W1);y=A(:,3);x=A(:,4,6);fort=1:Tt1=(t-1)*N+1;t2=t*N;wx(t1:t2,:)=W*x(t1:t2,:);endxconstant=ones(N*T,1);nobs K=size(x);info.lflag=0;% required f

22、or exact resultsinfo.model=1;info.fe=0;% Do not print intercept and fixed effects; use info.fe=1to turn onresults=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(r

23、esults,W,spat_model); panel_effects_sdm(results,vnames,W);183、时点固定效 应（Time period fixed effects）T=30;N=46;W=normw(W1);y=A(:,3);x=A(:,4,6);fort=1:Tt1=(t-1)*N+1;t2=t*N;wx(t1:t2,:)=W*x(t1:t2,:);endxconstant=ones(N*T,1);nobs K=size(x);info.lflag=0;% required for exact resultsinfo.model=2;info.fe=0;% Do

24、not print intercept and fixed effects; use info.fe=1to 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

25、;direct_indirect_effects_estimates(results,W,spat_model);19panel_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);fort=1:Tt1=(t-1)*N+1;t2=t*N;wx(t1:t2,:)=W*x(t1:t2,:);endxconstant=ones(N*T,1);nobs K=size(x);info.bc=0;info.lflag=0;%

26、required for exact resultsinfo.model=3;info.fe=0;% Do not print intercept and fixed effects; use info.fe=1to turn onresults=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_e

27、stimates(results,W,spat_model);20panel_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:KRafg(k,K+k)=1;% R(1,3)=0 and R(2,4)=0;endWald_spatial_lag=(Rafg*btemp)*inv(Rafg*

28、varcov*Rafg)*Rafg*btempprob_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:KR(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)+b

29、temp(4);endRafg=zeros(K,2*K+2);21fork=1:KRafg(k,k)=btemp(2*K+1);% k changed in 1, 7/12/2010Rafg(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);endWald_spatial_error=R*inv(Rafg*varcov*Rafg)*R prob_spatial_error=1

30、-chis_cdf (Wald_spatial_error,K)LR test spatial lagresultssar=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 errorresultssem=sem_panel_FE(y,x,W,T,info);LR_spatial_error=-2*(resultssem.lik-results.lik)prob_spatial

31、_error=1-chis_cdf (LR_spatial_error,K)5、空间随机效应与时点固定效 应模型22T=30;N=46;W=normw(W1);y=A(:,3);x=A(:,4,6);fort=1:Tt1=(t-1)*N+1;t2=t*N;wx(t1:t2,:)=W*x(t1:t2,:);endxconstant=ones(N*T,1);nobs K=size(x);ywith,xwith,meanny,meannx,meanty,meantx=demean(y,x wx,N,T,2); % 2=time dummiesinfo.model=1;results=sar_pane

32、l_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 lagbtemp=results.parm(1:2*K+2);23varcov=results.cov(1:2*K+2,1:2*K+2);Rafg=zeros(K,2*K+2);fork=1:KRafg(k,K+k)=1;% R(1,3)=0 an

33、d R(2,4)=0;endWald_spatial_lag=(Rafg*btemp)*inv(Rafg*varcov*Rafg)*Rafg*btempprob_spatial_lag= 1-chis_cdf (Wald_spatial_lag, K)wald test spatial errorR=zeros(K,1);fork=1:KR(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);

34、endRafg=zeros(K,2*K+2);fork=1:KRafg(k,k)=btemp(2*K+1);% k changed in 1, 7/12/2010Rafg(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);end24Wald_spatial_error=R*inv(Rafg*varcov*Rafg)*Rprob_spatial_error= 1-chis_cd

35、f (Wald_spatial_error,K)LR test spatial lagresultssar=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 errorresultssem=sem_panel_RE(ywith,xwith(:,1:K),W,T,info);LR_spatial_error=-2*(resultssem.lik-re

36、sults.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);fort=1:T25t1=(t-1)*N+1;t2=t*N;wx(t1:t2,:)=W*x(t1:t2,:);endxconstant=ones(N*T,1);nobs K=size(x);info.lflag=0;info.model=0;info.fe=0;results=sem_panel_FE(y,xconst

37、ant 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);fo

38、rt=1:Tt1=(t-1)*N+1;t2=t*N;wx(t1:t2,:)=W*x(t1:t2,:);26endxconstant=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_estimat

39、es(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);fort=1:Tt1=(t-1)*N+1;t2=t*N;wx(t1:t2,:)=W*x(t1:t2,:);endxconstant=ones(N*T,1);27nobs K=size(x);info.lflag=0;% required for exact resultsinfo.model=2;info.fe=0;%

40、 Do not print intercept and fixed effects; use info.fe=1to turn onresults=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

41、、双固定效应（Spatial and time period fixed effects）T=30;N=46;W=normw(W1);y=A(:,3);x=A(:,4,6);fort=1:Tt1=(t-1)*N+1;t2=t*N;wx(t1:t2,:)=W*x(t1:t2,:);endxconstant=ones(N*T,1);nobs K=size(x);28info.lflag=0;% required for exact resultsinfo.model=3;info.fe=0;% Do not print intercept and fixed effects; use info.f

42、e=1to turn onresults=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);fort=1:Tt1=(t-1)*N+1;t2=t*N;wx(t1:t2,:)=W*x(t1:t2,:);end29xconstant=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=strvca