两区制空间杜宾模型matlab代码

function results = sarregime_panel(y,x,d,W,T,info) % PURPOSE: computes two-regimes spatial lag model for spatial panels % (N regions*T time periods) with spatial fixed effects (u) and/or % time period fixed effects (v) % y = p1*W*y + p2*W*y + X*b + u (optional) + v (optional) + e, % using a binary variable d identifying the regimes % Supply data sorted first by time and then by spatial units, so first region 1, % region 2, et cetera, in the first year, then region 1, region 2, et % cetera in the second year, and so on % sarregime_panel computes y and x in deviation of the spatial and/or time means % --------------------------------------------------- % USAGE: results = sarregime_panel(y,x,d,W,T,info) % where: y = dependent variable vector % x = independent variables matrix % d = 0-1 variable identying to which regime each observations belongs % W = spatial weights matrix (standardized) % T = number of points in time % info = an (optional) structure variable with input options: % info.model = 0 pooled model without fixed effects (default, x may contain an intercept) % = 1 spatial fixed effects (x may not contain an intercept) % = 2 time period fixed effects (x may not contain an intercept) % = 3 spatial and time period fixed effects (x may not contain an intercept) % --------------------------------------------------- % RETURNS: a structure % results.meth = 'sarreg' if infomodel=0 % = 'sarregsfe' if info.model=1 % = 'sarregtfe' if info.model=2 % = 'sarregstfe' if info.model=3 % results.beta = bhat % results.rho = rho 1 and rho 2(p1 and p2 above) % results.tstat = asymp t-stat (last entries are rho1 and rho2=spatial autoregressive coefficients) % results.dif = rho1-rho2 % results.tdif = t-value of rho1-rho2 % results.cov = asymptotic variance-covariance matrix of the parameters b(eta) and rhos % results.yhat = yhat = [inv(I-p1*D*W-p2(I-D)*W2)]*[x*b+fixed effects] (according to prediction formula) % results.resid = y-p1*D*W*y-p2*(I-D)*W*y-x*b % results.sige = residuals'*residuals/nobs % results.rsqr = rsquared % results.corr2 = goodness-of-fit between actual and fitted values % results.sfe = spatial fixed effects (if info.model=1 or 3) % results.tfe = time period fixed effects (if info.model=2 or 3) % results.tsfe = t-values spatial fixed effects (if info.model=1 or 3) % results.ttfe = t-values time period fixed effects (if info.model=2 or 3) % results.con = intercept % results.con = t-value intercept % results.lik = log likelihood % results.nobs = # of observations % results.nvar = # of explanatory variables in x % results.tnvar = nvar + d*W*y + (I-D)*W*y + # fixed effects % results.iter = # of iterations taken % results.rmax = 1/max eigenvalue of W % results.rmin = 1/min eigenvalue of W % results.time = total time taken % -------------------------------------------------- % NOTE: Fixed effects and their t-values are calculated as the deviation % from the mean intercept % --------------------------------------------------- % % written by: J.Paul Elhorst 2/2007 % University of Groningen % Department of Economics % 9700AV Groningen % the Netherlands % j.p.elhorst@rug.nl % % REFERENCES: % Elhorst J.P., Fr�ret S. (2009) Evidence of political yardstick competition in France % using a two-regime spatial Durbin model with fixed effects. % Journal of Regional Science. Forthcoming. timet = clock; % start the clock for overall timing % if we have no options, invoke defaults if nargin == 5 info.model=0; fprintf(1,'default: pooled model without fixed effects \n'); end; fields = fieldnames(info); nf = length(fields); if nf > 0 for i=1:nf if strcmp(fields{i},'model') model = info.model; elseif strcmp(fields{i},'fe') fe = info.fe; end end end if model==0 results.meth='sarreg'; elseif model==1 results.meth='sarregsfe'; elseif model==2 results.meth='sarregtfe'; elseif model==3 results.meth='sarregstfe'; else error('sar_panel: wrong input number of info.model'); end % check size of user inputs for comformability [nobs nvar] = size(x); [N Ncol] = size(W); if N ~= Ncol error('sar: wrong size weight matrix W'); elseif N ~= nobs/T error('sar: wrong size weight matrix W or matrix x'); end; [nchk junk] = size(y); if nchk ~= nobs error('sar: wrong size vector y or matrix x'); end; lambda=eig(W); %eigenvalues rmin=min(lambda); rmax=max(lambda); % demeaning of the y and x variables, depending on (info.)model if (model==1 | model==3); meanny=zeros(N,1); meannx=zeros(N,nvar); for i=1:N ym=zeros(T,1); xm=zeros(T,nvar); for t=1:T ym(t)=y(i+(t-1)*N,1); xm(t,:)=x(i+(t-1)*N,:); end meanny(i)=mean(ym); meannx(i,:)=mean(xm); end clear ym xm; end % if statement if ( model==2 | model==3) meanty=zeros(T,1); meantx=zeros(T,nvar); for i=1:T t1=1+(i-1)*N;t2=i*N; ym=y([t1:t2],1); xm=x([t1:t2],:); meanty(i)=mean(ym); meantx(i,:)=mean(xm); end clear ym xm; end % if statement en=ones(T,1); et=ones(N,1); ent=ones(nobs,1); if model==1 ywith=y-kron(en,meanny); xwith=x-kron(en,meannx); elseif model==2 ywith=y-kron(meanty,et); xwith=x-kron(meantx,et); elseif model==3 ywith=y-kron(en,meanny)-kron(meanty,et)+kron(ent,mean(y)); xwith=x-kron(en,meannx)-kron(meantx,et)+kron(ent,mean(x)); else ywith=y; xwith=x; end % if statement % % Two regimes % index=d; Wy1=zeros(nobs,1); Wy2=zeros(nobs,1); wywith1=zeros(nobs,1); wywith2=zeros(nobs,1); for t=1:T t1=1+(t-1)*N;t2=t*N; W1=W; W2=W; for i=1:N if (index((t-1)*N+i,1)==0) W2(i,1:N)=zeros(1,N); else W1(i,1:N)=zeros(1,N); end end Wy1(t1:t2,1) = W1*y(t1:t2,1); Wy2(t1:t2,1) = W2*y(t1:t2,1); wywith1(t1:t2,1) = W1*ywith(t1:t2,1); wywith2(t1:t2,1) = W2*ywith(t1:t2,1); end % step 1) do regressions % step 2) maximize concentrated likelihood function; AI = xwith'*xwith; b0 = AI\(xwith'*ywith); bd1 = AI\(xwith'*wywith1); bd2 = AI\(xwith'*wywith2); e0 = ywith - xwith*b0; ed1 = wywith1 - xwith*bd1; ed2 = wywith2 - xwith*bd2; options = optimset('fminbnd'); rhos=[0.1;0.1]; [rhos,liktmp,exitflag,output]= fminsearch('f_sar2_panel',rhos,options,index,W,e0,ed1,ed2,N,T); if exitflag == 0 fprintf(1,'sar: convergence not obtained in %4d iterations \n',output.iterations); end; results.iter = output.iterations; % step 3) find b and sige maximum likelihood estimates rho1=rhos(1); rho2=rhos(2); beta = b0 - rho1*bd1- rho2*bd2; sige = (1/nobs)*(e0-rho1*ed1-rho2*ed2)'*(e0-rho1*ed1-rho2*ed2); results.rho = rhos; results.beta = beta; results.sige=sige; % step 4) find intercept and fixed effects if model==1 resh=y-rho1*Wy1-rho2*Wy2-x*beta; results.con=mean(resh); meanresn=zeros(N,1); for i=1:N rest=zeros(T,1); for t=1:T rest(t)=resh(i+(t-1)*N,1); end meanresn(i)=mean(rest); end clear rest; results.sfe=meanresn-kron(et,results.con); xhat=x*beta+kron(en,results.sfe)+kron(ent,results.con); tnvar=nvar+2+N; elseif model==2 resh=y-rho1*Wy1-rho2*Wy2-x*beta; results.con=mean(resh); meanrest=zeros(T,1); for t=1:T t1=1+(t-1)*N;t2=t*N; resn=resh(t1:t2,1); meanrest(t)=mean(resn); end clear resn; results.tfe=meanrest-kron(en,results.con); x