far_g.m

计量工具箱

function results = far_g(y,W,ndraw,nomit,prior)
% PURPOSE: Bayesian estimates for the 1st-order Spatial autoregressive model
%          y = rho*W*y + e,    e = N(0,sige*V), 
%          V = diag(v1,v2,...vn), r/vi = ID chi(r)/r, r = Gamma(m,k)
%          sige = gamma(nu,d0)    
%          rho = Uniform(rmin,rmax), or rho = beta(a1,a2); 
%----------------------------------------------------------------
% USAGE: result =  far_g(y,W,ndraw,nomit,prior)
% where: y = nobs x 1 independent variable vector (mean = 0)
%        W = nobs x nobs 1st-order contiguity matrix (standardized)
%       ndraw = # of draws
%       nomit = # of initial draws omitted for burn-in
%       prior = a structure variable for prior information input
%        prior.novi  = 1 turns off sampling for vi, producing homoscedastic model. default = 0            
%        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 (default = 1.01)
%        prior.a2    = (default = 1.01) see: 'help beta_prior'
%        prior.rval, = r prior hyperparameter, default=4
%        prior.m,    = informative Gamma(m,k) prior on r
%        prior.k,    = informative Gamma(m,k) prior on r
%        prior.eig   = 0 for default rmin = -1,rmax = +1, 1 for eigenvalue calculation of these
%        prior.rmin, = (optional) min value of rho to use in sampling
%        prior.rmax, = (optional) max value of rho to use in sampling
%        prior.lflag = 0 for full lndet computation (default = 1, fastest)
%                    = 1 for MC approximation (fast for very large problems)
%                    = 2 for Spline approximation (medium speed)
%        prior.dflag = 0 for numerical integration, 1 for Metropolis-Hastings (default = 0)
%        prior.order = order to use with info.lflag = 1 option (default = 50)
%        prior.iter  = iters to use with info.lflag = 1 option (default = 30)   
%        prior.lndet = a matrix returned by sar, far_g, sar_g, etc.
%                      containing log-determinant information to save time
%---------------------------------------------------------------
% RETURNS: a structure:
%          results.meth   = 'far_g'
%          results.pdraw  = rho draws (ndraw-nomit x 1)
%          results.sdraw  = sige draws (ndraw-nomit x 1)
%          results.vmean  = mean of vi draws (nobs x 1)
%          results.rdraw  = r-value draws (ndraw-nomit x 1)
%          results.nu     = prior nu-value for sige (if prior input)
%          results.d0     = prior d0-value for sige (if prior input)
%          results.a1     = a1 parameter for beta prior on rho from input, or default value
%          results.a2     = a2 parameter for beta prior on rho from input, or default value
%          results.r      = value of hyperparameter r (if input)
%          results.m      = m prior parameter (if input)
%          results.k      = k prior parameter (if input)    
%          results.nobs   = # of observations
%          results.ndraw  = # of draws
%          results.nomit  = # of initial draws omitted
%          results.y      = actual observations
%          results.yhat   = mean of posterior for y-predicted (nobs x 1)
%          results.time   = total time taken
%          results.time1  = time for log determinant calcluation
%          results.time2  = time for eigenvalue calculation   
%          results.time3  = time taken for sampling                 
%          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.dflag  = dflag value from input (or default value used)
%          results.iter   = 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.lndet = a matrix containing log-determinant information
%                          (for use in later function calls to save time)
%          results.acc   = acceptance rate for M-H sampling
%          results.mlike = log marginal likelihood (a vector ranging over
%                          rho values that can be integrated for model comparison)
%----------------------------------------------------------------
% NOTES: - use either improper prior.rval 
%          or informative Gamma prior.m, prior.k, not both of them
% - 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: (far_gd, far_gd2, far_gd3, demos) prt, plt
% --------------------------------------------------------------
% REFERENCES: James P. LeSage, `Bayesian Estimation of Spatial Autoregressive
%             Models',  International Regional Science Review, 1997 
%             Volume 20, number 1\&2, pp. 113-129.
% 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, last updated 8/2003
% Dept of Economics
% University of Toledo
% 2801 W. Bancroft St,
% Toledo, OH 43606
% [e-mail]

% 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.

timet = clock; % start the timer

if (nargin < 4 | nargin > 5)
    error('far_g: Wrong # of input arguments');
end;

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


 [n n2] = size(W);  
 if n ~= n2 % a non-square spatial weight matrix
 error('far_g: Wrong size 1st-order contiguity matrix');
 end; 
 [tst junk] = size(y);
 if tst ~= n % y-vector length doesn't match W-matrix
 error('far_g: Wrong size y vector on input');
 end;
 if junk ~= 1 % user didn't enter a column vector
 error('far_g: Wrong size y vector on input');
 end;

 
time1 = 0;
time2 = 0;
time3 = 0;

[nu,d0,rval,mm,kk,rho,sige,rmin,rmax,detval,ldetflag, ...
eflag,order,iter,novi_flag,cc,metflag,a1,a2] = far_parse(prior);

results.order = order;
results.iter = iter;



V = ones(n,1); in = ones(n,1); % initial value for V   
ys = y.*sqrt(V);
vi = in;
          
psave = zeros(ndraw-nomit,1);    % allocate storage for results
ssave = zeros(ndraw-nomit,1);
vmean = zeros(n,1);
yhat = zeros(n,1);
acc_rate = zeros(ndraw,1);


if mm~= 0                        % storage for draws on rvalue
rsave = zeros(ndraw-nomit,1);
end;

[rmin,rmax,time2] = far_eigs(eflag,W,rmin,rmax,n);

[detval,time1] = far_lndet(ldetflag,W,rmin,rmax,detval,order,iter);

iter = 1; 
time3 = clock; % start timing the sampler


% =====================================================
% The sampler starts here
% =====================================================

Wy = sparse(W)*y;
acc = 0;
cc = 0.1;

switch novi_flag    
case{0} % we do heteroscedastic model    
hwait = waitbar(0,'MCMC sampling far model ...');              

while (iter <= ndraw);        % start sampling;

% update sige;
nu1 = n + nu; 
Wys = Wy.*sqrt(V);
e = ys - rho*Wys;
d1 = d0 + e'*e;
chi = chis_rnd(1,nu1);
t2 = chi/d1;
sige = 1/t2;

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

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

     if metflag == 1
         % metropolis step to get rho update
          rhox = c_far(rho,y,sige,W,detval,in,a1,a2);
          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_far(rho2,y,sige,W,detval,in,a1,a2);
          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;
          rtmp(iter,1) = rho;
      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; % end of if metflag == 1


      if metflag == 0 % update rho using numerical integration
          e0 = ys;
          ed = Wys;
          epe0 = e0'*e0;
          eped = ed'*ed;
          epe0d = ed'*e0;
          
          rho = draw_rho(detval,epe0,eped,epe0d,n,1,rho,a1,a2);
      end;

    if (iter > nomit)
    ssave(iter-nomit,1) = sige;
    psave(iter-nomit,1) = rho;
    vmean = vmean + vi;
     if mm~= 0
     rsave(iter-nomit,1) = rval;
     end; 
    end; % end of if > nomit
    
iter = iter+1;
waitbar(iter/ndraw);

end; % end of sampling loop
close(hwait);

time3 = etime(clock,time3);


% compute posterior means and log marginal likelihood for return arguments
rho = mean(psave);
sigm = mean(ssave);
vmean = vmean/(ndraw-nomit);
V = in./vmean;

          ys = sqrt(V).*y;
          Wys = sqrt(V).*Wy;
          e0 = ys;
          ed = Wys;
          epe0 = e0'*e0;
          eped = ed'*ed;
          epe0d = ed'*e0;
e = (e0-rho*ed);
yhat = y - e;
sige = (1/n)*(e0-rho*ed)'*(e0-rho*ed); 
mlike = far_marginal(detval,e0,ed,epe0,eped,epe0d,n,1,a1,a2);


results.y = y;      
results.nobs = n;
results.nvar = 1;   
results.meth = 'far_g';
results.pdraw = psave;
results.sdraw = ssave;
results.vmean = vmean;
results.yhat = yhat;
results.resid = e;
results.tflag = 'plevel';
results.lflag = ldetflag;
results.dflag = metflag;
results.nobs  = n;
results.ndraw = ndraw;
results.nomit = nomit;
results.y = y;
results.nvar = 1;
results.mlike = mlike;
results.sige = sige;
results.rho = rho;
results.lndet = detval;
results.acc = acc_rate;
results.novi = novi_flag;

if mm~= 0
results.rdraw = rsave;
results.m     = mm;
results.k     = kk;
else
results.r     = rval;
results.rdraw = 0;
end;

results.time = etime(clock,timet);
results.time1 = time1;
results.time2 = time2;
results.time3 = time3;
results.lndet = detval;
results.rmax = rmax; 
results.rmin = rmin;

case{1} % we do homoscedastic model    
hwait = waitbar(0,'MCMC sampling far model ...');              

while (iter <= ndraw);        % start sampling;

% update sige;
nu1 = n + nu; 
e = y - rho*Wy;
d1 = d0 + e'*e;
chi = chis_rnd(1,nu1);
t2 = chi/d1;
sige = 1/t2;


     if metflag == 1
         % metropolis step to get rho update
          rhox = c_far(rho,y,sige,W,detval,in,a1,a2);
          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_far(rho2,y,sige,W,detval,in,a1,a2);
          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;
          rtmp(iter,1) = rho;
      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; % end of if metflag == 1