errors on extracting effects from dynamic spatial durbin model estimates

25 次查看（过去 30 天）
Henok Fasil Telila
Henok Fasil Telila 2019-9-24
回答： Henok Fasil Telila 2019-9-24
The following code is to estimate dynamic spatial durbin model and it works very well. EAVD, EVIC, and EVIP inside the code are dedicated to calculate the effects (one can check it around the middle of the code by this title "% Calculation of direct/indirect effects"). I tried to extract the direct, indirect and total effects from the result but I ain't successful. Any help please?
ylevelb=xlsread('SPATIAL_INDEX_RET_MALAB_NEW.xlsx','sheet 2'); %
x=xlsread('SPATIAL_VERY_FINAL_MATLAB1.xlsx','sheet 2');
WA=csvread('wm_exvol_NUOVO.csv');
T=226;
N=15;
ylevel=ylevelb;fprintf(1,'Y minus item 6 \n');
W=WA;fprintf(1,'W obv distance capitals \n');p=eig(W);eigmax=p(2);peig=1;
x3=x(:,[1:5,7:11]); %ERV
vnames3=strvcat('log(INDEX(t)+2)','log(INDEX(t-1)+2)','log(DEFAULT+2)','log(EXHRATE+2)','log(UNINFLA+2)','log(GDPGROWTH+2)','log(INTRATE+2)',...
'log(ylagerv+2)','log(slervcredT+2)','log(slervchecxT+2)','log(slervuninfT+2)','log(slervgdpT+2)','log(slervchintT+2)');
x=x3;vnames=vnames3;fprintf(1,'x var 1-7 + 9-12 \n');
%
% End choice of model
%
for t=1:T+1
t1=1+(t-1)*N;t2=t*N;
Wylevel(t1:t2,1)=W*ylevel(t1:t2,1);
end
y=ylevel(N+1:end)-ylevel(1:end-N);
%
% First the model without time-period fixed effects
%
info.lflag=0;
info.tl=1;
info.stl=1;
info.ted=1; %transformed approach
info.dyn=1;
info.model=1;
results=sar_panel_FE(ylevel(N+1:end),[ylevel(1:end-N) Wylevel(1:end-N) x],W,T,info);
prt_spnew(results,vnames,1);
results1=sar_jihai(ylevel,x,W,info);
results.beta=results1.theta1(1:end-2);
results.rho=results1.theta1(end-1);
results.tstat=results1.tstat(1:end-1);
results.sige=results1.theta1(end); %%%
results.lik = f2_sarpanel(results1.theta1,results1.yt,results1.zt,W,results.lndet,T); %%%
prt_spnew(results,vnames,1);
% Needed for F-test of time-period fixed effects
RRSS=results.sige*N*T;
%
% Now model with time-period fixed effect
%
info.lflag=0;
info.tl=1;
info.stl=1;
info.ted=1; %transformed approach
info.dyn=1;
info.model=3;
results=sar_panel_FE(ylevel(N+1:end),[ylevel(1:end-N) Wylevel(1:end-N) x],W,T,info);
%prt_spnew(results,vnames,1);
results1=sar_jihai_time(ylevel,x,W,info);
results.beta=results1.theta1(1:end-2);
results.rho=results1.theta1(end-1);
results.tstat=results1.tstat(1:end-1);
results.sige=results1.theta1(end); %%%
results.lik = f2_sarpanel(results1.theta1,results1.yt,results1.zt,W,results.lndet,T-1); %%%
prt_spnew(results,vnames,1);
varcov=results1.varcov;
R=zeros(1,1);
npar=length(results1.theta1);
tau = results1.theta1(1,1);
eta = results1.theta1(2,1);
rho = results1.theta1(npar-1,1);
R(1)=tau+rho+eta-1;
Rafg=zeros(1,npar);
Rafg(1,1)=1;Rafg(1,2)=1;Rafg(1,npar-1)=1;
sumpar=tau+rho+eta;
Wald=R(1)'*inv(Rafg(1,:)*varcov*Rafg(1,:)')*R(1);
F1=1-chis_cdf (Wald,1);
[sumpar Wald F1]
% F-test of time-period fixed effects
URSS=results.sige*N*T;
[Kjunck K1]=size([ylevel(1:end-N) Wylevel(1:end-N) x Wylevel(N+1:end)]);
F0=((RRSS-URSS)/(T-1))/(URSS/((N-1)*(T-1)-K1))
kansfo = fdis_prb(F0, T-1, (N-1)*(T-1)-K1)
%
% Spatial first-differences
%
sigma=(eye(N)-W)*(eye(N)-W)';
[V D]=eig(sigma);
transf=D(1+peig:end,1+peig:end)^(-0.5)*V(:,1+peig:end)'*(eye(N)-W);
Wstar=D(1+peig:end,1+peig:end)^(-0.5)*V(:,1+peig:end)'*W*V(:,1+peig:end)*D(1+peig:end,1+peig:end)^(0.5);
info.lflag=0;
for t=1:T+1
t1=1+(t-1)*N;t2=t*N;
ts1=1+(t-1)*(N-peig);ts2=t*(N-peig);
ystar(ts1:ts2,1)=transf*ylevel(t1:t2,1);
Wystar(ts1:ts2,1)=Wstar*ystar(ts1:ts2,1);
end
for t=1:T
t1=1+(t-1)*N;t2=t*N;
ts1=1+(t-1)*(N-peig);ts2=t*(N-peig);
xstar(ts1:ts2,:)=transf*x(t1:t2,:);
end
info.lflag=0;
info.tl=1;
info.stl=1;
info.dyn=1;
info.model=1;
N1=N-peig;
results=sar_panel_FE(ystar(N1+1:end),[ystar(1:end-N1) Wystar(1:end-N1) xstar],Wstar,T,info);
prt_spnew(results,vnames,1);
results1=sar_jihai(ystar,xstar,Wstar,info);
results.beta=results1.theta1(1:end-2);
results.rho=results1.theta1(end-1);
results.tstat=results1.tstat(1:end-1);
results.sige=results1.theta1(end); %%%
results.lik = f2_sarpanel(results1.theta1,results1.yt,results1.zt,Wstar,results.lndet,T-1); %%%
prt_spnew(results,vnames,1);
tau = results1.theta1(1,1);
eta = results1.theta1(2,1);
rho = results1.theta1(npar-1,1);
varcov=results1.varcov;
btemp=results1.theta1;
R=zeros(1,1);
R(1)=tau+rho+eta-1;
Rafg=zeros(1,npar);
Rafg(1,1)=1;Rafg(1,2)=1;Rafg(1,npar-1)=1;
sumpar=tau+rho+eta;
check=tau+eigmax*(rho+eta)-1; % this should be negative
Wald=R(1)'*inv(Rafg(1,:)*varcov*Rafg(1,:)')*R(1);
F1=1-chis_cdf (Wald,1);
[sumpar check Wald F1]
%
% Calculation of direct/indirect effects
%
NSIM=1000;
simresults=zeros(npar-2,NSIM);
simdir=zeros(npar-3,NSIM);
simind=zeros(npar-3,NSIM);
simtot=zeros(npar-3,NSIM);
for sim=1:NSIM
parms = chol(real(varcov)+0.001)'*randn(size(btemp)) + btemp;
rhosim = parms(npar-1,1);
betasim = parms(3:npar-2,1);
etasim=parms(1,1);
tausim = parms(2,1);
simresults(:,sim)=[tausim;betasim;rhosim];
SC=inv(eye(N)-rhosim*W)*((tausim-1)*eye(N)+(rhosim+etasim)*W);
p=1;
EAVD(p,1)=sum(diag(SC))/N; % average direct effect
EAVI(p,1)=sum(sum(SC,2)-diag(SC))/N; % average indirect effect
EAVC(p,1)=sum(sum(SC,1)'-diag(SC))/N; % average indirect effect
simdir(p,sim)=EAVD(p,1);
simind(p,sim)=EAVI(p,1);
simtot(p,sim)=EAVD(p,1)+EAVI(p,1);
for p=2:npar-3
C=zeros(N,N);
for i=1:N
for j=1:N
if (i==j) C(i,j)=betasim(p-1);
end
end
end
S=inv(eye(N)-rhosim*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
simdir(p,sim)=EAVD(p,1);
simind(p,sim)=EAVI(p,1);
simtot(p,sim)=EAVD(p,1)+EAVI(p,1);
end
end
[mean(simresults,2) mean(simresults,2)./std(simresults,0,2)]
fprintf(1,'First effect is convergence effect \n');
[mean(simdir,2) mean(simdir,2)./std(simdir,0,2) mean(simind,2) mean(simind,2)./std(simind,0,2)...
mean(simtot,2) mean(simtot,2)./std(simtot,0,2)]
%
% Impose restriction
%
theta2=btemp-varcov*Rafg'*inv(Rafg(1,:)*varcov*Rafg(1,:)')*R(1);
varcov2=varcov-varcov*Rafg'*inv(Rafg(1,:)*varcov*Rafg(1,:)')*Rafg*varcov;
results2=results;
results2.meth='sar_jihai_restricted';
results2.theta1=theta2;
results2.tstat1=theta2./(sqrt(abs(diag(varcov2))));
results2.sige=theta2(end); %%%
results2.lik = f2_sarpanel(results2.theta1,results1.yt,results1.zt,Wstar,results.lndet,T); %%%
help=theta2(end-1)*kron(speye(T),Wstar)*results1.yt;
residr=results1.yt-help-results1.zt*theta2(1:end-2);
yme=results1.yt-mean(results1.yt);
rsqr2=yme'*yme;
rsqr1 = residr'*residr;
results2.rsqr=1.0-rsqr1/rsqr2; %rsquared
res1=yme;
res2=((speye((N1)*T))-theta2(end-1)*kron(speye(T),Wstar))\(results1.zt*theta2(1:end-2))-mean(results1.yt);
rsq1=res1'*res2;
rsq2=res1'*res1;
rsq3=res2'*res2;
results2.corr2=rsq1^2/(rsq2*rsq3); %corr2
prt_sardynamic(results2,vnames,1);
tau = results2.theta1(1,1);
eta = results2.theta1(2,1);
rho = results2.theta1(npar-1,1);
%
% Calculation of direct/indirect effects
%
clear SC;
NSIM=1000;
simresults=zeros(npar-2,NSIM);
simdir=zeros(npar-3,NSIM);
simind=zeros(npar-3,NSIM);
simtot=zeros(npar-3,NSIM);
for sim=1:NSIM
parms = chol(real(varcov2)+0.001)'*randn(size(theta2)) + theta2;
rhosim = parms(npar-1,1);
betasim = parms(3:npar-2,1);
tausim = parms(2,1);
simresults(:,sim)=[tausim;betasim;rhosim];
SC=(tausim-1)*inv(eye(N)-rhosim*W)*(eye(N)-W);
p=1;
EAVD(p,1)=sum(diag(SC))/N; % average direct effect
EAVI(p,1)=sum(sum(SC,2)-diag(SC))/N; % average indirect effect
EAVC(p,1)=sum(sum(SC,1)'-diag(SC))/N; % average indirect effect
simdir(p,sim)=EAVD(p,1);
simind(p,sim)=EAVI(p,1);
simtot(p,sim)=EAVD(p,1)+EAVI(p,1);
for p=2:npar-3
C=zeros(N,N);
for i=1:N
for j=1:N
if (i==j) C(i,j)=betasim(p-1);
end
end
end
S=inv(eye(N)-rhosim*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
simdir(p,sim)=EAVD(p,1);
simind(p,sim)=EAVI(p,1);
simtot(p,sim)=EAVD(p,1)+EAVI(p,1);
end
end
[mean(simresults,2) mean(simresults,2)./std(simresults,0,2)]
fprintf(1,'First effect is convergence effect \n');
[mean(simdir,2) mean(simdir,2)./std(simdir,0,2) mean(simind,2) mean(simind,2)./std(simind,0,2)...
mean(simtot,2) mean(simtot,2)./std(simtot,0,2)]
clear simresults simdir simind simtot;
%end

请先登录，再回答此问题。

采纳的回答

Henok Fasil Telila
Henok Fasil Telila 2019-9-24
I did fix as follows. Use the code below and extract accordingly.
function direct_indirect_effects_estimates(results,W,spat_model)
N=results.N;
parm=results.parm;
cflag=results.cflag;
if (spat_model==0)
[junk nvar]=size(results.xwith);
NSIM=1000;
if (cflag==1) nvarc=nvar-1; else nvarc=nvar; end
simresults=zeros(nvarc+1,NSIM);
simdir=zeros(nvarc,NSIM);
simind=zeros(nvarc,NSIM);
simtot=zeros(nvarc,NSIM);
for sim=1:NSIM
parms = chol(results.cov)'*randn(size(parm)) + parm;
rhosim = parms(nvar+1,1);
if (results.cflag==1) betasim=parms(2:nvar,1);
else betasim=parms(1:nvar,1);
end
simresults(:,sim)=[betasim;rhosim];
for p=1:nvarc
C=zeros(N,N);
for i=1:N
for j=1:N
if (i==j) C(i,j)=betasim(p);
else C(i,j)=0;
end
end
end
S=inv(eye(N)-rhosim*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
simdir(p,sim)=EAVD(p,1);
simind(p,sim)=EAVI(p,1);
simtot(p,sim)=EAVD(p,1)+EAVI(p,1);
end
end
fprintf(1,' direct t-stat indirect t-stat total t-stat \n');
[mean(simdir,2) mean(simdir,2)./std(simdir,0,2) mean(simind,2) mean(simind,2)./std(simind,0,2)...
mean(simtot,2) mean(simtot,2)./std(simtot,0,2)]
elseif (spat_model==1)
[junk nvartot]=size(results.xwith);
if (cflag==1) nvar=(nvartot-1)/2; else nvar=nvartot/2; end
if (cflag==1) nvartotc=nvartot-1; else nvartotc=nvartot; end
NSIM=1000;
simresults=zeros(nvartotc+1,NSIM);
simdir=zeros(nvar,NSIM);
simind=zeros(nvar,NSIM);
simtot=zeros(nvar,NSIM);
for sim=1:NSIM
parms = chol(results.cov)'*randn(size(parm)) + parm;
rhosim = parms(nvartot+1,1);
if (results.cflag==1) betasim=parms(2:nvartot,1);
else betasim=parms(1:nvartot,1);
end
simresults(:,sim)=[betasim;rhosim];
for p=1:nvar
C=zeros(N,N);
for i=1:N
for j=1:N
if (i==j) C(i,j)=betasim(p);
else C(i,j)=betasim(nvar+p)*W(i,j);
end
end
end
S=inv(eye(N)-rhosim*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
simdir(p,sim)=EAVD(p,1);
simind(p,sim)=EAVI(p,1);
simtot(p,sim)=EAVD(p,1)+EAVI(p,1);
end
end
fprintf(1,' direct t-stat indirect t-stat total t-stat \n');
[mean(simdir,2) mean(simdir,2)./std(simdir,0,2) mean(simind,2) mean(simind,2)./std(simind,0,2)...
mean(simtot,2) mean(simtot,2)./std(simtot,0,2)]
else
error('wrong input number of spat_model');
end
%at the end run,
direct_indirect_effects_estimates(results,W,spat_model) ;
%results are stored regression values
%W is the weight matrix
%spat_model, is 1 if SDM and 0 if SAR

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Get Started with Model-Based Calibration Toolbox 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by