sarp_g.m

计量工具箱

function results = sarp_g(y,x,W,ndraw,nomit,prior)
% PURPOSE: Bayesian estimates of the spatial autoregressive probit model
%          y = rho*W*y + XB + e, e = N(0,I), 
%          B = N(c,T), 
%          rho = Uniform(rmin,rmax) 
%          y = binary, 0,1 variable
%-------------------------------------------------------------
% USAGE: results = sarp_g(y,x,W,ndraw,nomit,prior)
% where: y = dependent variable vector (nobs x 1)
%        x = independent variables matrix (nobs x nvar)
%        W = 1st order contiguity matrix (standardized, row-sums = 1)
%    ndraw = # of draws
%    nomit = # of initial draws omitted for burn-in            
%    prior = a structure variable with:
%            prior.novi  = 1 turns off sampling for vi, producing homoscedastic model  
%            prior.rval  = r prior hyperparameter, default = 4
%            prior.eig   = 0 for computing eigenvalues of W-matrix
%                          (defaults to 1, uses rmin = -1, rmax = 1)
%            prior.m     = informative Gamma(m,k) prior on r
%            prior.k     = (default: not used)
%            prior.nu    = informative Gamma(nu,d0) prior on sige
%            prior.d0    = default: nu=0,d0=0 (diffuse prior)
%            prior.a1    = parameter for beta(a1,a2) prior on rho see: 'help beta_prior'
%            prior.a2    = (default = 1.0, a uniform prior on rmin,rmax) 
%            prior.beta  = prior means for beta,   c above (default 0)
%            priov.bcov  = prior beta covariance , T above (default 1e+12)
%            prior.rmin  = (optional) min rho used in sampling (default = -1)
%            prior.rmax  = (optional) max rho used in sampling (default = 1)  
%            prior.lflag = 0 for full lndet computation (default = 1, fastest)
%                        = 1 for MC approx (fast for large problems)
%                        = 2 for Spline approx (medium speed)
%            prior.dflag = 0 for numerical integration, 1 for Metropolis-Hastings (default = 0)
%            prior.order = order to use with prior.lflag = 1 option (default = 50)
%            prior.iter  = iters to use with prior.lflag = 1 option (default = 30) 
%            prior.lndet = a matrix returned by sar, sar_g, sarp_g, etc.
%                          containing log-determinant information to save time
%-------------------------------------------------------------
% RETURNS:  a structure:
%          results.meth   = 'sarp_g'
%          results.bdraw  = bhat draws (ndraw-nomit x nvar)
%          results.pdraw  = rho  draws (ndraw-nomit x 1)
%          results.vmean  = mean of vi draws (nobs x 1) 
%          results.rdraw  = r draws (ndraw-nomit x 1) (if m,k input)
%          results.bmean  = b prior means, prior.beta from input
%          results.bstd   = b prior std deviations sqrt(diag(prior.bcov))
%          results.novi   = 1 for prior.novi = 1, 0 for prior.rval input
%          results.r      = value of hyperparameter r (if input)
%          results.nobs   = # of observations
%          results.nvar   = # of variables in x-matrix
%          results.ndraw  = # of draws
%          results.nomit  = # of initial draws omitted
%          results.nu     = nu prior parameter
%          results.d0     = d0 prior parameter
%          results.y      = y-vector from input (nobs x 1)
%          results.zip    = # of zero y-values
%          results.rsqr   = psuedo R-squared 
%          results.sige   = posterior mean of sige
%          results.yhat   = mean of posterior predicted (nobs x 1)
%          results.time1  = time for eigenvalue calculation
%          results.time2  = time for log determinant calcluation
%          results.time3  = time for sampling
%          results.time   = total time taken  
%          results.rmax   = 1/max eigenvalue of W (or rmax if input)
%          results.rmin   = 1/min eigenvalue of W (or rmin if input)          
%          results.tflag  = 'plevel' (default) for printing p-levels
%                         = 'tstat' for printing bogus t-statistics 
%          results.lflag  = lflag from input
%          results.iter   = prior.iter option from input
%          results.order  = prior.order option from input
%          results.limit  = matrix of [rho lower95,logdet approx, upper95] 
%                           intervals for the case of lflag = 1
%          results.acc    = an ndraw x 1 vector of acceptance rates for M-H sampling
%          results.lndet  = a matrix containing log-determinant information
%                          (for use in later function calls to save time)
%          results.mlike = log marginal likelihood (a vector ranging over
%                          rho values that can be integrated for model comparison)
% --------------------------------------------------------------
% NOTES: 
% - if you use lflag = 1 or 2, prior.rmin will be set = -1 
%                              prior.rmax will be set = 1
% - for n < 1000 you should use lflag = 0 to get exact results  
% --------------------------------------------------------------
% SEE ALSO: (sarp_gd, sarp_gd2 demos), prt
% --------------------------------------------------------------
% REFERENCES:  James P. LeSage, "Bayesian Estimation of Limited Dependent
%             variable Spatial Autoregressive Models", 
%             Geographical Analysis, 2000, Vol. 32, pp. 19-35.
% 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.
%----------------------------------------------------------------

% written by:
% James P. LeSage, 4/2002
% Dept of Economics
% University of Toledo
% 2801 W. Bancroft St,
% Toledo, OH 43606
% [e-mail]


% NOTE: some 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.

timet = clock;

% error checking on inputs
[n junk] = size(y);
[n1 k] = size(x);
[n3 n4] = size(W);
yin = y;
time1 = 0;
time2 = 0;
time3 = 0;

if n1 ~= n
error('sarp_g: x-matrix contains wrong # of observations');
elseif n3 ~= n4
error('sarp_g: W matrix is not square');
elseif n3~= n
error('sarp_g: W matrix is not the same size at y,x');
end;

if nargin == 5
    prior.lflag = 1;
end;

[nu,d0,rval,mm,kk,rho,sige,rmin,rmax,detval,ldetflag,eflag,order,iter,c,T,prior_beta,cc,metflag,novi_flag,a1,a2] = sar_parse(prior,k);

% error checking on prior information inputs
[checkk,junk] = size(c);
if checkk ~= k
error('sarp_g: prior means are wrong');
elseif junk ~= 1
error('sarp_g: prior means are wrong');
end;

[checkk junk] = size(T);
if checkk ~= k
error('sarp_g: prior bcov is wrong');
elseif junk ~= k
error('sarp_g: prior bcov is wrong');
end;

results.y = y;      
results.nobs = n;
results.nvar = k;   
results.order = order;
results.iter = iter;

timet = clock; % start the timer

[rmin,rmax,time1] = sar_eigs(eflag,W,rmin,rmax,n);

[detval,time2] = sar_lndet(ldetflag,W,rmin,rmax,detval,order,iter);

% storage for draws
          bsave = zeros(ndraw-nomit,k);
          psave = zeros(ndraw-nomit,1);
          ssave = zeros(ndraw-nomit,1);
          acc_rate = zeros(ndraw,1);
          ymean = zeros(n,1);
          yhat = zeros(n,1);
          yprob = zeros(n,1);

% ====== initializations
% compute this stuff once to save time
TI = inv(T);
TIc = TI*c;
iter = 1;

In = speye(n);
in = ones(n,1);
Wy = sparse(W)*y;
% find an index of values = 0
zipv = find(yin == 0);
zipo = find(yin == 1);
nzip = length(zipv);
sige = 1;
acc = 0;

switch novi_flag
    
case{0} % we do heteroscedastic model    
    
vmean = zeros(n,1);
if mm~= 0
    rsave = zeros(ndraw-nomit,1);
end;
W2diag = spdiags(W'*W,0);
V = ones(n,1);
vi = V;

hwait = waitbar(0,'sarp\_g: MCMC sampling ...');
t0 = clock;                  
iter = 1;
acc = 0;

          while (iter <= ndraw); % start sampling;
                  
          % update beta   
          xs = matmul(x,sqrt(V));
          ys = sqrt(V).*y;
          Wys = sqrt(V).*Wy;
          AI = inv(xs'*xs + sige*TI);		  
          yss = ys - rho*Wys;          
          b = xs'*yss + sige*TIc;
          b0 = AI*b;
          bhat = norm_rnd(sige*AI) + b0;  
          xb = xs*bhat;
          
          % update sige
          nu1 = n + 2*nu; 
          e = (yss - xb);
          d1 = 2*d0 + e'*e;
          chi = chis_rnd(1,nu1);
          sige = d1/chi;

          % update vi
          ev = y - rho*Wy - x*bhat; 
          chiv = chis_rnd(n,rval+1);   
          vi = ((ev.*ev/sige) + in*rval)./chiv;
          V = in./vi; 

          % update z-values
          mu = (In - rho*W)\(xb);
          ymu = y - mu;
           dsig = ones(n,1) + rho*rho*W2diag;
           yvar = ones(n,1)./dsig;
           A =  (1/sige)*(speye(n)-rho*W)*ymu; % a vector
           B =  (speye(n)-rho*W)'*A;  % a vector
           Cy = ymu - yvar.*B ;
           ym = mu + Cy;
          
           ind = find(yin == 0);
	       y(ind,1) = normrt_rnd(ym(ind,1),yvar(ind,1),0);
            
	       ind = find(yin == 1);
	       y(ind,1) = normlt_rnd(ym(ind,1),yvar(ind,1),0);
         
          % reformulate Wy
          Wy = sparse(W)*y;
         

              
          % update rval
          if mm ~= 0           
          rval = gamm_rnd(1,1,mm,kk);  
          end;

         if metflag == 1
         % metropolis step to get rho update
          rhox = c_sar(rho,ys,xb,sige,W,detval);
          accept = 0; 
          rho2 = rho + cc*randn(1,1); 
          while accept == 0
           if ((rho2 > rmin) & (rho2 < rmax)); 
           accept = 1;  
           else
           rho2 = rho + cc*randn(1,1);
           end; 
          end;
          rhoy = c_sar(rho2,ys,xb,sige,W,detval);
          ru = unif_rnd(1,0,1);
          if ((rhoy - rhox) > exp(1)),
          p = 1;
          else, 
          ratio = exp(rhoy-rhox); 
          p = min(1,ratio);
          end;
              if (ru < p)
              rho = rho2;
              acc = acc + 1;
              end;
      acc_rate(iter,1) = acc/iter;
      % update cc based on std of rho draws
       if acc_rate(iter,1) < 0.4
       cc = cc/1.1;
       end;
       if acc_rate(iter,1) > 0.6
       cc = cc*1.1;
       end;
    end;

         if metflag == 0
      % when metflag == 0,
      % we use numerical integration to perform rho-draw
          b0 = (xs'*xs)\(xs'*ys);
          bd = (xs'*xs)\(xs'*Wys);
          e0 = ys - xs*b0;
          ed = Wys - xs*bd;
          epe0 = e0'*e0;
          eped = ed'*ed;
          epe0d = ed'*e0;
          logdetx = log(det(xs'*xs));
          rho = draw_rho(detval,epe0,eped,epe0d,n,k,rho,a1,a2,logdetx);
      end;



               
    if iter > nomit % if we are past burn-in, save the draws
    bsave(iter-nomit,1:k) = bhat';
    psave(iter-nomit,1) = rho;
    ssave(iter-nomit,1) = sige;
    ymean = ymean + y;
    vmean = vmean + vi;
        if mm~= 0
        rsave(iter-nomit,1) = rval;
        end;         
    end;
    
iter = iter + 1; 
waitbar(iter/ndraw);         
end; % end of sampling loop
close(hwait);

% compute posterior means and log marginal likelihood for return arguments
ymean = ymean /(ndraw-nomit);
bmean = mean(bsave);
beta = bmean';
rho = mean(psave);
vmean = vmean/(ndraw-nomit);
results.vmean = vmean;
V = in./vmean;
yhat = (speye(n) - rho*W)\(x*beta);
yprob = stdn_cdf(yhat);

          Wy = W*ymean;
          xs = matmul(x,sqrt(V));
          ys = sqrt(V).*ymean;
          Wys = sqrt(V).*Wy;
          AI = inv(xs'*xs);
          b0 = (xs'*xs)\(xs'*ys);
          bd = (xs'*xs)\(xs'*Wys);
          e0 = ys - xs*b0;
          ed = Wys - xs*bd;
          epe0 = e0'*e0;
          eped = ed'*ed;
          epe0d = ed'*e0;

sige = (1/(n-results.nvar))*(e0-rho*ed)'*(e0-rho*ed); 
logdetx = log(det(xs'*xs));
[nobs,nvar] = size(xs);
mlike = sar_marginal(detval,e0,ed,epe0,eped,epe0d,nobs,nvar,logdetx,a1,a2);

% compute psuedo R-squared
e = ymean-yhat;
sigu = (e'*e);
sige = sigu/(nobs-nvar);
ym = ymean - mean(ymean);
rsqr1 = sigu;
rsqr2 = ym'*ym;
rsqr = 1.0 - rsqr1/rsqr2; % psuedo r-squared

case{1} % homoscedastic model


hwait = waitbar(0,'sarp\_g: MCMC sampling ...');
t0 = clock;                  
iter = 1;
acc = 0;
W2diag = spdiags(W'*W,0);

          while (iter <= ndraw); % start sampling;
                  
          % update beta   
          AI = inv(x'*x + sige*TI);
          ys = y - rho*Wy;          
          b = x'*ys + sige*TIc;
          bm = AI*b;
          bhat = norm_rnd(sige*AI) + bm;  
          xb = x*bhat;
          
          % update sige
          nu1 = n + 2*nu; 
          e = (ys - xb);
          d1 = 2*d0 + e'*e;
          chi = chis_rnd(1,nu1);
          sige = d1/chi;

          % update z-values
          mu = (In - rho*W)\(xb);
          ymu = y - mu;
           dsig = ones(n,1) + rho*rho*W2diag;
           yvar = ones(n,1)./dsig;
           A =  (1/sige)*(speye(n)-rho*W)*ymu; % a vector
           B =  (speye(n)-rho*W)'*A;  % a vector
           Cy = ymu - yvar.*B ;
           ym = mu + Cy;
          
           ind = find(yin == 0);
	       y(ind,1) = normrt_rnd(ym(ind,1),yvar(ind,1),0);
            
	       ind = find(yin == 1);
	       y(ind,1) = normlt_rnd(ym(ind,1),yvar(ind,1),0);
         
          % reformulate Wy
          Wy = sparse(W)*y;
         
              

         if metflag == 1
         % metropolis step to get rho update
          rhox = c_sar(rho,y,xb,sige,W,detval);
          accept = 0; 
          rho2 = rho + cc*randn(1,1); 
          while accept == 0
           if ((rho2 > rmin) & (rho2 < rmax)); 
           accept = 1;  
           else
           rho2 = rho + cc*randn(1,1);
           end; 
          end;
          rhoy = c_sar(rho2,y,xb,sige,W,detval);
          ru = unif_rnd(1,0,1);
          if ((rhoy - rhox) > exp(1)),
          p = 1;
          else, 
          ratio = exp(rhoy-rhox); 
          p = min(1,ratio);
          end;