程序：
% dimensions of the problem
A=xlsread('C:\Users\lenovo\Documents\MATLAB\Zandberg\cigarette.xls');
W1=xlsread('C:\Users\lenovo\Documents\MATLAB\Zandberg\Spat-Sym-US.xls');
T=29; % 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+1
    t1=(t-1)*N+1;t2=t*N;
    wx(t1:t2,:)=W*x(t1:t2,:);
    Wy(t1:t2,1)=W*y(t1:t2,1);
end
[nobs K]=size(x);
xconstant=ones(nobs,1);
info.lflag=0;
info.tl=1;
info.stl=1;
info.ted=1;
info.dyn=1;
info.model=3;
info.fe=0;
info.bc=0;
results=sar_panel_FE(y(N+1:end),[y(1:end-N) Wy(1:end-N) x(N+1:end,:) wx(N+1:end,:)],W,T,info);
vnames=strvcat('logcit','timelag logcit','spacetimelag logcit','logp','logy','W*logp','W*logy');
results1=sar_jihai_time(y(1:nobs),[x(N+1:nobs,:) wx(N+1:nobs,:)],W,info);
results.beta=results1.theta1(1:end-2);
results.rho=results1.theta1(end-1);
results.tstat=results1.tstat(1:end-1);
prt_sp(results,vnames,1);
btemp=results1.theta1;
varcov=results1.varcov;
%
% Direct and indirect effects estimates
% st=short term, lt=long term
% c=(long term) convergence effect of dependent variable
%
% Direct and indirect effects estimates
% st=short term, lt=long term
% c=(long term) convergence effect of dependent variable
%
NSIM=1000;
px=2;
[npar dummy]=size(btemp);
simresults=zeros(npar-1,NSIM);
simdirst=zeros(px,NSIM);
simindst=zeros(px,NSIM);
simtotst=zeros(px,NSIM);
simdirlt=zeros(px,NSIM);
simindlt=zeros(px,NSIM);
simtotlt=zeros(px,NSIM);
simdirc=zeros(1,NSIM);
simindc=zeros(1,NSIM);
simtotc=zeros(1,NSIM);
for sim=1:NSIM
    parms = chol(varcov)'*randn(size(btemp)) + btemp;
    deltasim = parms(npar-1,1); % coef WY(t)
    betasim = parms(3:npar-2,1);
    tausim = parms(1,1); % Coef Y(t-1)
    etasim = parms(2,1); % Coef WY(t-1)
    simresults(:,sim)=[tausim;etasim;betasim;deltasim];
    SS=(eye(N)-deltasim*W)\eye(N);
    SC=SS*((tausim-1)*eye(N)+(deltasim+etasim)*W);
    simdirc(1,sim)=sum(diag(SC))/N; % average direct effect
    simindc(1,sim)=sum(sum(SC,2)-diag(SC))/N; % average indirect effect
    simtotc(1,sim)=simdirc(1,sim)+simindc(1,sim);
    for p=1:px
        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(p+2)*W(i,j);
                end
            end
        end
        SC=SS*C;
        simdirst(p,sim)=sum(diag(SC))/N; % average direct effect
        simindst(p,sim)=sum(sum(SC,2)-diag(SC))/N; % average indirect effect
        simtotst(p,sim)=simdirst(p,sim)+simindst(p,sim);
        SC=((1-tausim)*eye(N)-(deltasim+etasim)*W)\C;
        simdirlt(p,sim)=sum(diag(SC))/N; % average direct effect
        simindlt(p,sim)=sum(sum(SC,2)-diag(SC))/N; % average indirect effect
        simtotlt(p,sim)=simdirlt(p,sim)+simindlt(p,sim);        
    end
end
fprintf(1,'Convergence effect \n');
[mean(simdirc,2) mean(simdirc,2)./std(simdirc,0,2) mean(simindc,2) mean(simindc,2)./std(simindc,0,2)...
    mean(simtotc,2) mean(simtotc,2)./std(simtotc,0,2)]
fprintf(1,'Short term effects \n');
[mean(simdirst,2) mean(simdirst,2)./std(simdirst,0,2) mean(simindst,2) mean(simindst,2)./std(simindst,0,2)...
    mean(simtotst,2) mean(simtotst,2)./std(simtotst,0,2)]
fprintf(1,'Long term effects \n');
[mean(simdirlt,2) mean(simdirlt,2)./std(simdirlt,0,2) mean(simindlt,2) mean(simindlt,2)./std(simindlt,0,2)...
    mean(simtotlt,2) mean(simtotlt,2)./std(simtotlt,0,2)]
clear simresults simdirc simindc simtotc simdirst simindst simtotst simdirlt simindlt simtotlt;

运行结果：

Pooled model with spatially lagged dependent variable, spatial and time period fixed effects
Dependent Variable = logcit              
R-squared          =    0.9772   
corr-squared       =    0.8615   
sigma^2            =    0.0011
Nobs,Nvar,#FE      =   1334,     7,    80  
log-likelihood     =        2623.2622
# of iterations    =      1   
min and max rho    =   -1.3924,   1.0000
total time in secs =    0.4510
time for optimiz   =    0.1110
time for lndet     =    0.0130
time for eigs      =    0.1290
time for t-stats   =    0.0080
No lndet approximation used
***************************************************************
Variable                 Coefficient  Asymptot t-stat    z-probability
timelag logcit              0.864568        65.041225         0.000000
spacetimelag logcit        -0.015347        -0.287520         0.773714
logp                       -0.265979       -13.192203         0.000000
logy                        0.099796         4.165387         0.000031
W*logp                      0.170413         3.664979         0.000247
W*logy                     -0.022369        -0.868904         0.384900
W*dep.var.                  0.076402         2.025253         0.042841

??? Reference to non-existent field 'varcov'.

Error in ==> demo at 33
varcov=results1.varcov;