dynamic_spatial_panel_动态空间计量模型代码_Able_

function results = sar_jihai_time(y,x,W,info) % PURPOSE: computes spatial autoregressive model estimates % Ynt = p*Wn*Ynt +gamma*Yn,t-1+rho*WnYn,t-1+ X*b+c +at+ e, using sparse matrix algorithms % --------------------------------------------------- % USAGE: results = sar(y,x,W,info) % where: y = dependent variable vector % x = explanatory variables matrix % W = standardized contiguity matrix % c=individual fixed effects % at=time dummy effects % info = an (optional) structure variable with input options: % info.rmin = (optional) minimum value of rho to use in search (default = -1) % info.rmax = (optional) maximum value of rho to use in search (default = +1) % info.eig = 0 for default rmin = -1,rmax = +1, 1 for eigenvalue calculation of these % info.convg = (optional) convergence criterion (default = 1e-8) % info.maxit = (optional) maximum # of iterations (default = 500) % info.lflag = 0 for full lndet computation (default = 1, fastest) % = 1 for MC lndet approximation (fast for very large problems) % = 2 for Spline lndet approximation (medium speed) % info.order = order to use with info.lflag = 1 option (default = 50) % info.iter = iterations to use with info.lflag = 1 option (default = 30) % info.lndet = a matrix returned by sar, sar_g, sarp_g, etc. % containing log-determinant information to save time % --------------------------------------------------- % RETURNS: a structure % results.meth = 'sar_jihai_time' % results.theta = estimates vector before bias correction with the % order gammma, rho, b, p, sigma % results.theta1 = estimates vector after bias correction with the order gammma, rho, b, p, sigma % results.tstat = asymp t-stat (last entry is rho) % results.bstd = std of betas (nvar x 1) vector % results.pstd = std of rho % results.yhat = yhat (nobs x 1) vector % results.resid = residuals (nobs x 1) vector % results.sige = sige % results.rsqr = rsquared % results.rbar = rbar-squared % results.lik = log likelihood % results.nobs = # of observations % results.nvar = # of explanatory variables in x % results.y = y data vector % results.iter = # of iterations taken % results.rmax = 1/max eigenvalue of W (or rmax if input) % results.rmin = 1/min eigenvalue of W (or rmin if input) % results.lflag = lflag from input % results.liter = info.iter option from input % results.order = info.order option from input % results.limit = matrix of [rho lower95,logdet approx, upper95] intervals % for the case of lflag = 1 % results.time1 = time for log determinant calcluation % results.time2 = time for eigenvalue calculation % results.time3 = time for hessian or information matrix calculation % results.time4 = time for optimization % results.time = total time taken % results.lndet = a matrix containing log-determinant information % (for use in later function calls to save time) % -------------------------------------------------- % NOTES: if you use lflag = 1 or 2, info.rmin will be set = -1 % info.rmax will be set = 1 % For n < 1000 you should use lflag = 0 to get exact results % -------------------------------------------------- % REFERENCES: Anselin (1988), pages 180-182. % For lndet information see: Ronald Barry and R. Kelley Pace, % "A Monte Carlo Estimator of the Log Determinant of Large Sparse Matrices", % Linear Algebra and its Applications", Volume 289, Number 1-3, 1999, pp. 41-54. % and: R. Kelley Pace and Ronald P. Barry % "Simulating Mixed Regressive Spatially autoregressive Estimators", % Computational Statistics, 1998, Vol. 13, pp. 397-418. % --------------------------------------------------- % extended by Jihai Yu, 05/23/2007, from the original code % written by: % James P. LeSage, 1/2000 % Dept of Economics % University of Toledo % 2801 W. Bancroft St, % Toledo, OH 43606 % jlesage@spatial.econometrics.com % NOTE: much of the speed for large problems comes from: % the use of methods pioneered by Pace and Barry. % R. Kelley Pace was kind enough to provide functions % lndetmc, and lndetint from his spatial statistics toolbox % for which I'm very grateful. fields = fieldnames(info); nf = length(fields); for i=1:nf, if strcmp(fields{i},'tl') tl = info.tl; % tar lag index elseif strcmp(fields{i},'stl') stl = info.stl; % star lag index elseif strcmp(fields{i},'ted') ted = info.ted; % ted=1 transformation approach, =2, direct approachs end end [junk n]=size(W); [tempsize junk]=size(y); t=tempsize/n-1; if ted==1 %use the transformation approach Jn=speye(n)-1/n*ones(n,1)*ones(1,n); [Fnn junk]=eig(Jn); F=Fnn(:,2:n); Fy=kron(speye(t+1),F'); Fx=kron(speye(t),F'); y2=Fy*y; if isempty(x) == 0 x2=Fx*x; else x2=[]; end W2=F'*W*F; n=n-1; info2=struct('n',n,'t',t,'rmin',0,'rmax',1,'lflag',0,'tl',tl,'stl',stl); res=sar_jihai(y2,x2,W2,info2); results.meth = 'sar_jihai_time'; results.yt = res.yt; results.nobs = res.nobs; results.rmax = res.rmax; results.rmin = res.rmin; results.order = res.order; results.miter = res.miter; results.time1 = res.time1; results.time2 = res.time2; results.time3 = res.time3; results.time4 = res.time4; results.lndet = res.lndet; results.lik = res.lik; results.kz = res.kz; results.varcov = res.varcov; results.bias1=res.bias1; results.SIG=res.SIG; results.SIGi=res.SIGi; results.OMG=res.OMG; results.SIG1=res.SIG1; results.SIGi1=res.SIGi1; results.OMG1=res.OMG1; results.theta=res.theta; results.theta1=res.theta1; results.tstat=res.tstat; results.tstat1=res.tstat1; results.std=res.std; results.std1=res.std1; results.rsqr=res.rsqr; %Added J.P.Elhorst, 20-4-2010 results.corr2=res.corr2; %Added J.P.Elhorst, 20-4-2010 results.zt=res.zt; %Added J.P.Elhorst, 20-4-2010 results.respaul=res.respaul; %Added J.P.Elhorst, 12-1-2016 results.respaul1=res.respaul1; %Added J.P.Elhorst, 12-1-2016 results.lik1=res.lik1; %Added J.P.Elhorst, 31-3-2016 results.resid=res.resid; %Added J.P.Elhorst, 31-3-2016 else if ted==0 error('Wrong Info input,use sar_jihai instead as there is no time dummy'); else %ted==2, use the direct approach Wnt=kron(speye(t),W); nt=n*t; yt=y(n+1:n+nt); ytl=y(1:nt); Q=kron(eye(t)-1/t*ones(t,1)*ones(1,t),eye(n)-1/n*ones(n,1)*ones(1,n)); yt=Q*yt; ytl=Q*ytl; ysl=Wnt*yt; ystl=Wnt*ytl; xt=x; if isempty(x) == 0 xt=Q*xt; else xt=[]; end if stl + tl == 2 zt=[ytl ystl xt]; elseif stl + tl == 1 if stl == 1, zt=[ystl xt]; else zt=[ytl xt]; end elseif stl + tl == 0 error('Wrong Info input,Our model has dynamic term anyway'); else error('Double-Check stl & tl # in Info structure '); end [junk kz]=size(zt); [junk kx]=size(xt); time1 = 0; time2 = 0; time3 = 0; time4 = 0; timet = clock; % start the clock for overall timing