errors on extracting effects from dynamic spatial durbin model estimates

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

Answer 1

[Henok Fasil Telila]

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