目前用matlab做空间计量还处于讳莫如深的境况，很多会操作的人还利用它来赚钱（诸如开培训班等），而不是把相关代码的收集和运算公开，让更多人受益。我觉得这都没有必要，我们做科学研究不是窝里争，再者很多代码是国外学者编的，所以没必要藏着掖着，而是让科学流行起来。在这之前，我公开了R和stata的相关命令和操作，都很受欢迎。基于我的理解和掌握，我决定贡献并适当讲解一下matlab代码，本次分享的代码是：空间面板数据，可能还不是很周全，抛砖引玉吧！如有不懂的地方，请到Regional Study群交流（234081931）。
- A=xlsread('D:\matlab\toolbox\SpatEcon_ren\Files_SLX_paper\cigarette.xls');
- W1=xlsread('D:\matlab\toolbox\SpatEcon_ren\Files_SLX_paper\Spat-Sym-US.xls');
- %
- % SAR, SEM, SAC, SDM, SDEM and GNS models. Every time this file is run
- %
- T=30; % number of time periods
- N=46; % number of regions
- % row-normalize W
- W=normw(W1); % function of LeSage
- y=A(:,[3]); % column number in the data matrix that corresponds to the dependent variable
- x=A(:,[4,6]); % column numbers in the data matrix that correspond to the independent variables
- 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); %size返回的是行和列
- et=ones(T,1);
- en=ones(N,1);
- %
- % All models include spatial and time period fixed effects
- %
- % ----------------------------------------------------------------------------------------
- % No spatial interaction effects
- 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);
- 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);
- loglik=-nobs/2*log(2*pi*sige)-1/(2*sige)*results.resid'*results.resid
- % ----------------------------------------------------------------------------------------
- % Spatial lag model
- 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
- % New routines to calculate effects estimates
- results=sar_panel_FE(y,x,W,T,info);
- vnames=strvcat('logcit','logp','logy');
- % Print out coefficient estimates
- prt_sp(results,vnames,1);
- % Print out effects estimates
- spat_model=0;
- direct_indirect_effects_estimates(results,W,spat_model);
- % panel_effects_sar(results,vnames,W);
- % ----------------------------------------------------------------------------------------
- % Spatial error model
- 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
- % New routines to calculate effects estimates
- results=sem_panel_FE(y,x,W,T,info);
- % Print out coefficient estimates
- prt_sp(results,vnames,1);
- % ----------------------------------------------------------------------------------------
- % Model with exogenous interaction effects
- model=3;
- [ywith,xwith,meanny,meannx,meanty,meantx]=demean(y,[x wx],N,T,model);
- results=ols(ywith,xwith);
- vnames=strvcat('logcit','logp','logy','W*logp','W*logy');
- prt_reg(results,vnames);
- intercept=mean(y)-mean([x wx])*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 wx]*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);
- loglik=-nobs/2*log(2*pi*sige)-1/(2*sige)*results.resid'*results.resid
- % ----------------------------------------------------------------------------------------
- % Spatial Durbin model
- 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
- info.bc=1;
- % New routines to calculate effects estimates
- results=sar_panel_FE(y,[x wx],W,T,info);
- % Print out coefficient estimates
- prt_sp(results,vnames,1);
- % Print out effects estimates
- spat_model=1;
- direct_indirect_effects_estimates(results,W,spat_model);
- %panel_effects_sdm(results,vnames,W);
- % ----------------------------------------------------------------------------------------
- % Spatial Durbin error model
- 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
- info.bc=1;
- % New routines to calculate effects estimates
- results=sem_panel_FE(y,[x wx],W,T,info);
- % Print out coefficient estimates
- prt_sp(results,vnames,1);
- % ----------------------------------------------------------------------------------------
- % SAC / SARAR / Clifford-Ord/ Kelejian-Prucha model
- 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
- info.bc=1;
- results=sac_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);
- % ----------------------------------------------------------------------------------------
- % Full model
- info.model=3;
- info.fe=0; % Do not print intercept and fixed effects; use info.fe=1 to turn on
- info.bc=1;
- results=sac_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);
Pual Elhost的主页
http://www.regroningen.nl/elhorst/software.shtml
James Lesage 个人主页
http://www.spatial-econometrics.com/
kelley pace个人网站
http://www.spatial-statistics.com/
Donald J. Lacombe个人主页
http://community.wvu.edu/~djlacombe/
Journal of Statistical Software杂志首页
http://www.jstatsoft.org/ |
熟悉论坛请点击新手指南
京公网安备 11010802022788号
论坛法律顾问：王进律师
知识产权保护声明
免责及隐私声明
