mess_g2.m

计量工具箱

function results = mess_g2(y,x,options,ndraw,nomit,prior,start)
% PURPOSE: Bayesian estimates of the matrix exponential spatial model (mess)
% % [samples values of rho to produce a posterior distribution]
% S*y = X*b + e,           with xflag == 0, or:
% S*y = [i X D*X]*b + e,   with xflag == 1
%          e = N(0,sige*In), 
%          S = e^alpha*D
%          B = N(c,T), 
%          1/sige = Gamma(nu,d0), 
%          alpha = Uniform(amin,amax) or alpha = N(a,B) 
% D = a weight matrix constructed from neighbors N_i using:
% D = sum rho^i N_i / sum rho^i, i=1,...,#neighbors 
% a grid of rho values from options.rmin to options.rmax
% a single value for # neighbors, option.neigh is used in this function
%-------------------------------------------------------------
% USAGE: results = mess_g2(y,x,options,ndraw,nomit,prior,start)
% where: y = dependent variable vector (nobs x 1)
%        x = independent variables matrix (nobs x nvar)
%  options = a structure variable with:
%  options.latt  = lattitude coordinates (nx1 vector)
%  options.long  = longitude coordinates (nx1 vector)
%  options.neigh = # of neighbors to use in constructing D
%  options.xflag = 0 for S*y = X*b + e,         model (default)
%                = 1 for S*y = [i X D*X]*b + e, model
%  options.rmin  = minimum value of rho to use in discounting 
%                  (0 < rho < 1), (default 0.5)
%  options.rmax  = maximum value of rho to use in discounting (default 1)
%  options.nflag = 0 for neighbors using 1st and 2nd order Delauney (default)
%                = 1 for neighbors using 3rd and 4th order Delauney
%                  for large nobs, and large # neighbors, used nflag = 1
%  options.q     = # of terms to use in the matrix exponential
%                  expansion (default = 7)
%    ndraw = # of draws
%    nomit = # of initial draws omitted for burn-in            
%    prior = a structure variable with:
%            prior.beta  = prior means for beta,   c above (default 0)
%                          (nvar x 1 vector, D*X terms have diffuse prior =0)
%            priov.bcov  = prior beta covariance , T above (default 1e+12)
%                          (nvar x nvar matrix, D*X terms have diffuse prior)
%            prior.alpha = prior mean for alpha    a above (default uniform)
%            prior.rcov  = prior alpha variance,   B above 
%            prior.nu    = informative Gamma(nu,d0) prior on sige
%            prior.d0    = default: nu=0,d0=0 (diffuse prior)
%            prior.m,    = informative Gamma(m,k) prior on r
%            prior.k,    = informative Gamma(m,k) prior on r
%            prior.amin  = (optional) min alpha used in sampling 
%                          (default = max like alpha - 3*std(alpha))
%            prior.amax  = (optional) max alpha used in sampling 
%                          (default = max like alpha + 3*std(alpha))
%            prior.lflag = 0 (default for marginal likelihood calculations)
%                        = 1 for no marginal likelihood (faster)
%    start = (optional) structure containing starting values: 
%            defaults: beta=ones(k,1),sige=1,rho=0.5, V=ones(n,1)
%            start.b   = beta starting values (nvar x 1)
%            start.a   = alpha starting value (scalar)
%            start.sig = sige starting value  (scalar)
%-------------------------------------------------------------
% RETURNS:  a structure:
%          results.meth   = 'mess_g2'
%          results.bdraw  = bhat draws (ndraw-nomit x nvar)
%          results.bmean  = mean of bhat draws
%          results.bstd   = std of bhat draws
%          results.adraw  = alpha draws (ndraw-nomit x 1)
%          results.amean  = mean of alpha draws
%          results.astd   = std of alpha draws
%          results.sdraw  = sige draws (ndraw-nomit x 1)
%          results.smean  = mean of sige draws
%          results.lmean  = marginal likelihood based on mean of draws
%          results.rmean  = posterior mean of rho
%          results.rstd   = posterior std of rho
%          results.rdraw  = draws for rho
%          results.bprior = b prior means, prior.beta from input
%          results.bpstd  = b prior std deviations sqrt(diag(prior.bcov))
%          results.nobs   = # of observations
%          results.nvar   = # of variables in x-matrix (plus D*X matrix)
%          results.ndraw  = # of draws
%          results.nomit  = # of initial draws omitted
%          results.y      = y-vector from input (nobs x 1)
%          results.yhat   = mean of posterior predicted (nobs x 1)
%          results.nu     = nu prior parameter
%          results.d0     = d0 prior parameter
%          results.stime  = time for sampling
%          results.time   = total time taken  
%          results.ntime  = time taken for mesh over rho and alpha values
%          results.accept = acceptance rate 
%          results.amax   = amax: max like alpha + 3*std(alpha) (or user input)
%          results.amin   = amin: max like alpha - 3*std(alpha) (or user input)      
%          results.rmax   = rmax: default (or user input)
%          results.rmin   = rmin: default (or user input)        
%          results.tflag  = 'plevel' (default) for printing p-levels
%                         = 'tstat' for printing bogus t-statistics 
%          results.palpha   = prior for alpha (from input)
%          results.acov   = prior variance for alpha (from input)
%          results.pflag  = 1, if a normal(a,B) prior for alpha, 0 otherwise
%          results.xflag = model flag from input
%          results.neigh = # of terms in flexible D-matrix specification
%                         (from input or default)
%          results.q     = q value from input (or default)

% --------------------------------------------------------------
% NOTES: if the model includes a constant term
% it should be entered as the first column in the x-matrix
% that is input to the function
% 1) mess_g1 produces a posterior distribution for # neighbors
% 2) mess_g2 produces a posterior distribution for the hyperparameter rho
% 3) mess_g3 produces posteriors for both rho and # of neighbors
% --------------------------------------------------------------
% SEE ALSO:  mess_g2d, mess_g, mess_g3, messv_g, prt, c_mess, mess
% --------------------------------------------------------------
% REFERENCES: LeSage and Pace (2000) "Bayesian Estimation of the
% Matrix Exponential Spatial Specification", unpublished manuscript
%----------------------------------------------------------------

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


timet = clock;

% error checking on inputs
[n junk] = size(y);
results.y = y;
[n1 k] = size(x);

if n1 ~= n
error('mess_g2: x-matrix contains wrong # of observations');
end;

% set defaults
q = 7;
xflag = 0;
rmin = 0.5;
rmax = 1.0;
neigh = 0;
aflag = 0;
llflag = 0;
nflag = 0;
pflag = 0; % flag for the presence or absent of a prior on alpha

mm = 0;    % set defaults
nu = 0;    % default diffuse prior for sige
d0 = 0;
sig0 = 1;         % default starting values for sige
astart = -1;      % default starting value for alpha
c = zeros(k,1);   % diffuse prior for beta
T = eye(k)*1e+12;
palpha = -1;
S = 1e+12;
lflag = 0; % default to do marginal likelihood calculations
if nargin == 7
    if ~isstruct(start)
        error('mess_g2: must supply starting values in a structure');
    end;
 % parse starting values entered by the user
 fields = fieldnames(start);
 nf = length(fields);
 for i=1:nf
    if strcmp(fields{i},'b')
        b0 = start.b; 
        [n1 n2] = size(b0); % error checking on user inputs
       if n1 ~= k
        error('mess_g2: starting beta values are wrong');
       elseif n2 ~= 1
        error('mess_g2: starting beta values are wrong');
       end;
    elseif strcmp(fields{i},'sig')
        sig0 = start.sig;
       [n1 n2] = size(sig0); % error checking on user inputs
       if n1 ~= 1
        error('mess_g2: starting sige value is wrong');
       elseif n2 ~= 1
        error('mess_g2: starting sige value is wrong');
       end;
    elseif strcmp(fields{i},'a')
        astart = start.a;
       [n1 n2] = size(astart); % error checking on user inputs
       if n1 ~= 1
        error('mess_g2: starting alpha value is wrong');
       elseif n2 ~= 1
        error('mess_g2: starting alpha value is wrong');
       end;
    end;
 end; % end of for loop
% parse options structure
    if ~isstruct(options)
        error('mess_g2: must supply option values in a structure');
    end;
 fields = fieldnames(options);
 nf = length(fields);
 for i=1:nf
    if strcmp(fields{i},'xflag')
       xflag = options.xflag;
    elseif strcmp(fields{i},'rmin')
        rmin = options.rmin;
    elseif strcmp(fields{i},'rmax')
        rmax = options.rmax;
    elseif strcmp(fields{i},'q')
       q = options.q;
     elseif strcmp(fields{i},'neigh')
       neigh = options.neigh;
     elseif strcmp(fields{i},'latt')
        latt = options.latt; llflag = llflag + 1;
     elseif strcmp(fields{i},'long')
        long = options.long; llflag = llflag + 1;
    elseif strcmp(fields{i},'nflag')
        nflag = options.nflag;
     end;
 end; % end of for loop

% parse prior structure variable inputs        
    if ~isstruct(prior)
    error('mess_g2: must supply the prior as a structure variable');
    end;

fields = fieldnames(prior);
nf = length(fields);

for i=1:nf
    if strcmp(fields{i},'beta')
        c = prior.beta;
    elseif strcmp(fields{i},'bcov')
        T = prior.bcov;
    elseif strcmp(fields{i},'alpha')
        palpha = prior.alpha; pflag = 1;
    elseif strcmp(fields{i},'acov')
        S = prior.acov;        
    elseif strcmp(fields{i},'nu')
        nu = prior.nu;
    elseif strcmp(fields{i},'d0')
        d0 = prior.d0;
    elseif strcmp(fields{i},'lflag')
       lflag = prior.lflag; 
    end;
end;


elseif nargin == 6   % we supply default starting values
 fields = fieldnames(prior);
 nf = length(fields);
 for i=1:nf
    if strcmp(fields{i},'beta')
        c = prior.beta;
    elseif strcmp(fields{i},'bcov')
        T = prior.bcov;
    elseif strcmp(fields{i},'alpha')
        palpha = prior.alpha; pflag = 1;
    elseif strcmp(fields{i},'acov')
        S = prior.acov;         
    elseif strcmp(fields{i},'nu')
        nu = prior.nu;
    elseif strcmp(fields{i},'d0')
        d0 = prior.d0;
    elseif strcmp(fields{i},'rval')
       rval = prior.rval; 
    elseif strcmp(fields{i},'lflag')
       lflag = prior.lflag; 
    end;
 end;
 
 % parse options
     if ~isstruct(options)
        error('mess_g2: must supply option values in a structure');
    end;
 fields = fieldnames(options);
 nf = length(fields);
 for i=1:nf
    if strcmp(fields{i},'xflag')
       xflag = options.xflag;
    elseif strcmp(fields{i},'rmin')
        rmin = options.rmin;
    elseif strcmp(fields{i},'rmax')
        rmax = options.rmax;
    elseif strcmp(fields{i},'q')
       q = options.q;
     elseif strcmp(fields{i},'neigh')
       neigh = options.neigh;
     elseif strcmp(fields{i},'latt')
        latt = options.latt; llflag = llflag + 1;
     elseif strcmp(fields{i},'long')
        long = options.long; llflag = llflag + 1;
    elseif strcmp(fields{i},'nflag')
        nflag = options.nflag;
     end;
 end; % end of for loop



elseif nargin == 5   % we supply all defaults
   % parse options structure
    if ~isstruct(options)
        error('mess_g2: must supply option values in a structure');
    end;
 fields = fieldnames(options);
 nf = length(fields);
 for i=1:nf
    if strcmp(fields{i},'xflag')
       xflag = options.xflag;
    elseif strcmp(fields{i},'rmin')
        rmin = options.rmin;
    elseif strcmp(fields{i},'rmax')
        rmax = options.rmax;
    elseif strcmp(fields{i},'q')
       q = options.q;
     elseif strcmp(fields{i},'neigh')
       neigh = options.neigh;
     elseif strcmp(fields{i},'latt')
        latt = options.latt; llflag = llflag + 1;
     elseif strcmp(fields{i},'long')
        long = options.long; llflag = llflag + 1;
    elseif strcmp(fields{i},'nflag')
        nflag = options.nflag;
     end;
 end; % end of for loop

else
error('Wrong # of arguments to mess_g2');
end;
      

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

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

[checkk junk] = size(palpha);
if checkk ~= 1
error('mess_g2: prior alpha is wrong');
elseif junk ~= 1
error('mess_g2: prior alpha is wrong');
end;

[checkk junk] = size(S);
if checkk ~= 1
error('mess_g2: prior acov is wrong');
elseif junk ~= 1
error('mess_g: prior acov is wrong');
end;

% make sure the user input latt, long or we really bomb
if llflag ~= 2;
error('mess_g2: no lattitude-longitude coordinates input');
end;

switch xflag % switch on x transformation

   
   case{0} % case where x variables are not transformed


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

% ========= do up front grid over rho, alpha values

% find maxl alpha using ? rho value ( 1 used here, the default);

results.rmin = rmin;
results.rmax = rmax;

rgrid = rmin:0.01:rmax;
nrho = length(rgrid);