用R语言呈现Elhost书中的表2.2

28关注
175
粉丝

已卖：437份资源

学科带头人

4%

还不是VIP/贵宾

-

0%

威望: 1 级
论坛币: 4428 个
通用积分: 61.8219
学术水平: 89 点
热心指数: 97 点
信用等级: 67 点
经验: 1609 点
帖子: 715
精华: 1
在线时间: 1963 小时
注册时间: 2011-7-21
最后登录: 2025-5-15

楼主

jianhui80

发表于 2016-12-19 21:55:55 |AI写论文

是否 +2 论坛币

k人参与回答

经管之家送您一份

应届毕业生专属福利!

求职就业群

赵安豆老师微信：zhaoandou666

经管之家联合CDA

送您一个全额奖学金名额~ !

立即领取

感谢您参与论坛问题回答

经管之家送您两个论坛币！

+2 论坛币

由Elhost编写，肖光恩译的《空间计量经济学：从横截面到面板》的matlab代码可以到Elhost个人主页上查找，下面由我给出R语言代码，今天给出的是表2.2的结果。需要申明的是，由于stargazer包现在只能识别errorsarlm和lagsarlm，故而在表没有显示SLX和GNS。

library(spdep)
library(stargazer)
data(oldcol)
lw <- nb2listw(COL.nb, style="W")
OLS <- lm(CRIME ~ INC + HOVAL , data = COL.OLD)
SAR <- lagsarlm(CRIME ~ INC + HOVAL, data=COL.OLD,nb2listw(COL.nb, style="W"), method="eigen", quiet=FALSE)
SEM <- errorsarlm(CRIME ~ INC + HOVAL, data=COL.OLD,lw, method="eigen", quiet=FALSE)
SLX <- lmSLX(CRIME ~ INC + HOVAL, data=COL.OLD, listw=lw)
SAC <- sacsarlm(CRIME ~ INC + HOVAL, data=COL.OLD,nb2listw(COL.nb, style="W"))
SDM <- lagsarlm(CRIME ~ INC + HOVAL, data=COL.OLD,nb2listw(COL.nb, style="W"), method="eigen", quiet=FALSE,type = "mixed")
SDEM <- errorsarlm(CRIME ~ INC + HOVAL, data=COL.OLD,lw, method="eigen", quiet=FALSE,etype = "emixed")
GNS <- sacsarlm(CRIME ~ INC + HOVAL, data=COL.OLD,nb2listw(COL.nb, style="W"),type = "sacmixed")
stargazer(OLS,SAR,SEM,SDM,SDEM, type = "text",model.names = F, column.labels = c("OLS","SAR","SEM","SDM","SDEM"), title=" comparison of the estimation results explaining the crime rate")

复制代码

注:（1）截距项是原来的100倍，标准误稍有出入；
（2）对数似然值差别很大，这个还需找到计算上的差别。

扫码加我拉你入群

请注明：姓名-公司-职位

以便审核进群资格，未注明则拒绝

分享0 收藏6 回帖

关键词：host HOS R语言 matlab代码 Library library 经济学 matlab 横截面空间

相关帖子

区域经济学是经济学中还未开放的花骨朵。

沙发

xixia333 发表于 2016-12-19 23:10:11

很赞，希望楼主继续分享，谢谢！

藤椅

wingate99 发表于 2016-12-20 14:41:47

感谢楼主！！！！

板凳

jianhui80

发表于 2017-1-10 17:32:48

% Anselin's Columbus neighborhood spatial data sample
%
%
load wmat.dat; load anselin.dat; % load the data;
anselin(:,1:3)=anselin(:,1:3)/100;
y = anselin(:,1);
W = wmat;
N=49;
W=zeros(N,N);
for p=1:232
W(wmat(p,1),wmat(p,2))=1;
end
W=normw(W);
lambda=eig(W);
nobs = length(y);
% construct an explanatory variables matrix
x = [ones(nobs,1) anselin(:,2) anselin(:,3)];
wx = W*anselin(:,2:3);
% set up variable names
vnames = strvcat('Crime','constant','income','house value');
% OLS model
results = ols(y,x); % compare to least-squares estimates
prt(results,vnames); % print results
loglikelihood = -(N/2)*log(2*pi) - (N/2)*log(results.si2) - (1/results.si2)*results.resid'*results.resid
LMsarsem_panel(results,W,y,x); % (Robust) LM tests
lm1=lmlag_panel(y,[ones(nobs,1) x],W);
prt_tests(lm1);
lm2=lmerror_panel(y,[ones(nobs,1) x],W);
prt_tests(lm2);
lm3=lmlag_robust_panel(y,[ones(nobs,1) x],W);
prt_tests(lm3);
lm4=lmerror_robust_panel(y,[ones(nobs,1) x],W);
prt_tests(lm4);
% Spatial lag model
result = sarpaul(y,x,W); % estimate spatial autoregressive model
prt(result,vnames); % print results
% Exact direct/indirect effects estimates
nvar=2;
beta=result.beta(2:3);
rho=result.rho;
for p=1:nvar
C=zeros(N,N);
for i=1:N
for j=1:N
if (i==j) C(i,j)=beta(p);
% else C(i,j)=beta(nvar+p)*W(i,j);
end
end
end
S=(eye(N)-rho*W)\C;
EAVD(p,1)=sum(diag(S))/N; % average direct effect
EAVI(p,1)=sum(sum(S,2)-diag(S))/N; % average indirect effect
EAVC(p,1)=sum(sum(S,1)'-diag(S))/N; % average indirect effect
EAVCtot(p,1)=EAVD(p,1)+EAVI(p,1);
end
fprintf(1,' direct indirect total \n');
[EAVD EAVI EAVCtot]
result.N=result.nobs;
result.cflag=1;
result.xwith=x;
spat_model=0;
direct_indirect_effects_estimates(result,W,spat_model);
%Spatial error model
result = sempaul(y,x,W); % estimate spatial error model
prt(result,vnames); % print results
% OLS model with exogenous interaction effects
vnames = strvcat('Crime','constant','income','house value','W*income','W*house value');
results = ols(y,[x wx]); % compare to least-squares estimates
prt(results,vnames); % print results
loglikelihood = -(N/2)*log(2*pi) - (N/2)*log(results.si2) - (1/results.si2)*results.resid'*results.resid
% Spatial Durbin model
vnames = strvcat('Crime','constant','income','house value','W*income','W*house value');
v1names = strvcat('Crime','constant','income','house value');
result = sarpaul(y,[x wx],W); % estimate spatial Durbin model
prt(result,vnames); % print results
% Exact direct/indirect effects estimates
nvar=2;
beta=result.beta(2:5);
rho=result.rho;
for p=1:nvar
C=zeros(N,N);
for i=1:N
for j=1:N
if (i==j) C(i,j)=beta(p);
else C(i,j)=beta(nvar+p)*W(i,j);
end
end
end
S=(eye(N)-rho*W)\C;
EAVD(p,1)=sum(diag(S))/N; % average direct effect
EAVI(p,1)=sum(sum(S,2)-diag(S))/N; % average indirect effect
EAVC(p,1)=sum(sum(S,1)'-diag(S))/N; % average indirect effect
EAVCtot(p,1)=EAVD(p,1)+EAVI(p,1);
end
fprintf(1,' direct indirect total \n');
[EAVD EAVI EAVCtot]
result.N=result.nobs;
result.cflag=1;
result.xwith=[x wx];
spat_model=1;
direct_indirect_effects_estimates(result,W,spat_model);
% Spatial Durbin error model
result = sempaul(y,[x wx],W); % estimate spatial error model
prt(result,vnames); % print results
% SAC / SARAR / Clifford-Ord/ Kelejian-Prucha model
result = sacpaul(y,x,W);
prt(result,v1names); % print results
% Exact direct/indirect effects estimates
nvar=2;
beta=result.beta(2:3);
rho=result.rho;
for p=1:nvar
C=zeros(N,N);
for i=1:N
for j=1:N
if (i==j) C(i,j)=beta(p);
% else C(i,j)=beta(nvar+p)*W(i,j);
end
end
end
S=(eye(N)-rho*W)\C;
EAVD(p,1)=sum(diag(S))/N; % average direct effect
EAVI(p,1)=sum(sum(S,2)-diag(S))/N; % average indirect effect
EAVC(p,1)=sum(sum(S,1)'-diag(S))/N; % average indirect effect
EAVCtot(p,1)=EAVD(p,1)+EAVI(p,1);
end
fprintf(1,' direct indirect total \n');
[EAVD EAVI EAVCtot]
result.N=result.nobs;
result.cflag=1;
result.xwith=x;
spat_model=0;
direct_indirect_effects_estimates(result,W,spat_model);
%%% Full model
result = sacpaul(y,[x wx],W);
prt(result,vnames); % print results
parm=[result.beta;result.rho;result.lam;result.sige];
% Exact direct/indirect effects estimates
nvar=2;
beta=result.beta(2:5);
rho=result.rho;
for p=1:nvar
C=zeros(N,N);
for i=1:N
for j=1:N
if (i==j) C(i,j)=beta(p);
else C(i,j)=beta(nvar+p)*W(i,j);
end
end
end
S=(eye(N)-rho*W)\C;
EAVD(p,1)=sum(diag(S))/N; % average direct effect
EAVI(p,1)=sum(sum(S,2)-diag(S))/N; % average indirect effect
EAVC(p,1)=sum(sum(S,1)'-diag(S))/N; % average indirect effect
EAVCtot(p,1)=EAVD(p,1)+EAVI(p,1);
end
fprintf(1,' direct indirect total \n');
[EAVD EAVI EAVCtot]
result.parm=parm;
result.N=result.nobs;
result.cflag=1;
result.xwith=[x wx];
spat_model=1;
direct_indirect_effects_estimates(result,W,spat_model);

复制代码

报纸

jianhui80

发表于 2017-1-10 17:36:04

上面的Matlab代码来自Elhost主页CRSEModels文件中的demo_crime_rates.m

地板

jianhui80

发表于 2017-1-10 17:54:34

我们可以在R代码的基础上得到OLS的AIC和对数似然值
logLik_lm.model
Log likelihood of the linear model for rho=0
AIC_lm.model
AIC of the linear model for rho=0

SAR$logLik_lm.model
SAR$AIC_lm.model

复制代码

7楼

jianhui80

发表于 2017-1-10 18:01:57

当然也可以直接用

logLik(OLS)
AIC(OLS)

复制代码

8楼

jianhui80

发表于 2017-1-10 18:38:45

尽管我们发现用R和Matlab得到的对数似然值完全不一样，然而，我们发现似然比却是一致的，比如书中35页中写到“OLS模型的对数似然值从13.776增加到17.075。对SLX和OLS模型进行似然比检验的统计量为6.598，其自由度为2，在5%的水平下临界值为6”。
LR.sarlm(SLX,OLS)

      Likelihood ratio for spatial linear models

data:
Likelihood ratio = 6.5981, df = 2, p-value = 0.03692
sample estimates:
Log likelihood of SLX Log likelihood of OLS
         -184.0782          -187.3772

9楼

w财经w

发表于 2018-1-11 20:58:08

jianhui80 发表于 2017-1-10 17:32

您好，楼主大神，在matlab上运行您给的这个程序，出现了
??? Reference to non-existent field 'si2'.

Error in ==> RENda at 25
loglikelihood = -(N/2)*log(2*pi) - (N/2)*log(results.si2) -
(1/results.si2)*results.resid'*results.resid;

求教怎么破？？？

10楼

w财经w

发表于 2018-1-12 08:58:59

jianhui80 发表于 2017-1-10 17:32

楼主，上个提问的问题自己解决了，问下从您第27行开始的稳健LM估计出现了以下问题：
??? Undefined function or method 'prt_tests' for input arguments of type 'struct'.

Error in ==> demo_crime_rates at 27
prt_tests(lm1);

求教如何解决，自己改了半天程序都不行啊

[区域经济学] 用R语言呈现Elhost书中的表2.2 [推广有奖]

经管之家送您一份

经管之家联合CDA

感谢您参与论坛问题回答

扫码加我拉你入群

相关帖子

浏览过的帖子

浏览过的版块

本版微信群

[区域经济学] 用R语言呈现Elhost书中的表2.2 [推广有奖]

经管之家送您一份

经管之家联合CDA

感谢您参与论坛问题回答

扫码加我 拉你入群

相关帖子

浏览过的帖子

浏览过的版块

本版微信群

扫码加我拉你入群