messv_g3.m

计量工具箱

          v(i,1) = alpha.^(i-1);
          end;
          W = (1./[1 cumprod(1:nq1)]);
          Sy = Ycap*diag(W)*v;
          Sys = Ycaps*diag(W)*v;

          % update beta   
          AI = inv(Xcaps'*Xcaps + sige*TI);
          b = Xcaps'*Sys + sige*TIc;
          b0 = AI*b;
          bhat = norm_rnd(sige*AI) + b0;  
         
          % update sige
          nu1 = n + 2*nu; 
          e = (Sys - Xcaps*bhat);
          d1 = 2*d0 + e'*e;
          chi = chis_rnd(1,nu1);
          sige = d1/chi;

          % update vi
          e = Sy - Xcap*bhat;
          chiv = chis_rnd(n,rval+1);   
          vi = ((e.*e./sige) + in*rval)./chiv;
          V = in./vi; 
  
          % update rval
          if mm ~= 0           
          rval = gamm_rnd(1,1,mm,kk);  
          end;


          % metropolis step to get alpha update
          if pflag == 0
          alphax = cn_mess3(alpha,y,Xcap,Ymat,bhat,sige,rho,neigh,rmat,mmat); 
          elseif pflag == 1
          alphax = cn_mess3(alpha,y,Xcap,Ymat,bhat,sige,rho,neigh,rmat,mmat,palpha,S); 
          end;
          accept = 0;
          alpha2 = alpha + cc*randn(1,1);
          while accept == 0
           if (alpha2 <= 0) 
           accept = 1;  
           else
           alpha2 = alpha + cc*randn(1,1);
           cntr = cntr+1; % counts accept rate
           end; 
          end; 
           if pflag == 0
           alphay = cn_mess3(alpha2,y,Xcap,Ymat,bhat,sige,rho,neigh,rmat,mmat);
           elseif pflag == 1
           alphay = cn_mess3(alpha2,y,Xcap,Ymat,bhat,sige,rho,neigh,rmat,mmat,palpha,S);
           end;
          ru = unif_rnd(1,0,1);
          if ((alphay - alphax) > exp(1)),
          p = 1;
          else,          
          ratio = exp(alphay-alphax);
          p = min(1,ratio);
          end;
              if (ru < p)
                 alpha = alpha2;
              end;

          rtmp(iter,1) = alpha;

          % update rho using metroplis-hastings step
          rhox = rn_mess3(rho,y,Xcap,Ymat,bhat,alpha,neigh,rmat,mmat); 
          rho2 = unif_rnd(1,rmin,rmax);
          rhoy = rn_mess3(rho2,y,Xcap,Ymat,bhat,alpha,neigh,rmat,mmat);
           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;
              end;
	  
          % update neigh using metroplis-hastings step
          neighx = nn_mess3(neigh,y,Xcap,Ymat,bhat,alpha,rho,rmat,mmat); 
          neigh2 = round(unif_rnd(1,mmin,mmax));

          neighy = nn_mess3(neigh2,y,Xcap,Ymat,bhat,alpha,rho,rmat,mmat);
           ru = unif_rnd(1,0,1);
          if ((neighy - neighx) > exp(1)),
          p = 1;
          else,          
          ratio = exp(neighy-neighx);
          p = min(1,ratio);
          end;
              if (ru < p)
                 neigh = neigh2;
              end;

	  
     % evaulate the likelihood using current draws
     if lflag == 0
                like = -(n/2)*log(2*pi*sige) - (e'*e)/(2*sige);
     end;
                       
     % update rval
     if mm ~= 0           
     rval = gamm_rnd(1,1,mm,kk);  
     end;
              
    if iter > nomit % if we are past burn-in, save the draws
    bsave(iter-nomit,:) = bhat';
    ssave(iter-nomit,1) = sige;
    asave(iter-nomit,1) = alpha; 
    rsave(iter-nomit,1) = rho;
    msave(iter-nomit,1) = neigh;
    if mm~= 0
    rdraw(iter-nomit,1) = rval;
    end;
    vmean = vmean + vi;

    if lflag == 0
       lsave = lsave + like;
    else
       lsave = lsave + 0;
    end;
    
    
    end;

       if iter == nomit % update cc based on initial draws
         tst = 2*std(rtmp(1:nomit,1));
         if tst > 0.05
         cc = tst;
         end;
    end;
                 

iter = iter + 1; 
waitbar(iter/ndraw);         
end; % end of sampling loop
close(hwait);

stime = etime(clock,t0);

% compute posterior means
if lflag == 0
lmean = lsave/(ndraw-nomit);
else 
   lmean = 0;
end;
amean = mean(asave);
bmean = mean(bsave);
astd = std(asave);
bstd = std(bsave);
smean = mean(ssave);
rmean = mean(rsave);
rstd = std(rsave);
mmean = mean(msave);
mstd = std(msave);
vmean = vmean/(ndraw-nomit);

% find acceptance rate
results.accept = 1 - cntr/(iter+cntr);

mround = round(mmean);

tmp = rmean.^(0:mround-1);
tmp = tmp/sum(tmp);

nnlist = nnlistall(:,1:mround);

% construct Sy, Sx based on posterior means
wy = y;
Y = y(:,ones(1,q));
for i=2:q;
wy = wy(nnlist)*tmp';
Y(:,i) = wy;
end;
Ys = matmul(Y,sqrt(vmean));
[junk nq] = size(Y);
nq1 = nq-1;
v = ones(nq,1);
for i=2:nq;
v(i,1) = amean.^(i-1);
end;
W = (1./[1 cumprod(1:nq1)]);
sy = Y*diag(W)*v;
sys = Ys*diag(W)*v;

% create  Sx based on posterior mean of rho
xmat = x;
for i=1:nk;
xi = xsub(:,i);
tmpp = xi(nnlist)*tmp';
xmat = [xmat tmpp];
end;
xmats = matmul(xmat,sqrt(vmean));

e = sy - xmat*bmean';
yhat = y - e;
sigu = e'*e;
ym = y - mean(y);
rsqr1 = sigu;
rsqr2 = ym'*ym;
rsqr = 1.0 - rsqr1/rsqr2; % r-squared
rsqr1 = rsqr1/(n-k);
rsqr2 = rsqr2/(n-1.0);
rbar = 1 - (rsqr1/rsqr2); % rbar-squared

time = etime(clock,timet);

results.meth  = 'messv_g3';
results.bdraw = bsave;
results.adraw = asave;
results.bmean = bmean';
results.bstd  = bstd';
results.amean = amean;
results.astd  = astd;
results.smean = smean;
results.sdraw = ssave;
results.lmean = lmean;
results.rmean = rmean;
results.rstd  = rstd;
results.rdraw = rsave;
results.mmean = mmean;
results.mstd  = mstd;
results.mdraw = msave;
results.vmean = vmean;
results.rvdraw = rdraw;
results.bprior = c;
results.bpstd  = sqrt(diag(Tnew));
results.nobs  = n;
results.nvar  = k;
results.ndraw = ndraw;
results.nomit = nomit;
results.time  = time;
results.stime = stime;
results.nu = nu;
results.d0 = d0;
if mm~= 0
results.m = m;
results.k = k;
else
results.rval = rval;
end;
results.tflag = 'plevel';
results.pflag = pflag;
results.palpha = palpha;
results.acov = S;
results.y = y;
results.yhat = yhat;
results.resid = e;
results.rsqr = rsqr;
results.rbar = rbar;
results.q     = q;
results.ntime = gtime;
results.nobs = n;
results.nvar = k;
results.xflag = xflag;
results.nflag = nflag;

otherwise
   
end; % end of switch


function cout = rn_mess3(rho,ys,xs,Ymat,beta,alpha,neigh,rmat,mmat);
% PURPOSE: evaluate the conditional distribution of alpha 
%          for the Bayesian mess_g3 model
% ---------------------------------------------------
%  USAGE: cout = c_mess3(alpha,y,x,Symat,beta,alpha,rho,neigh,amat,rmat,mmat)
%  where:  alpha  = matrix exponential alpha spatial parameter
%            y    = dependent variable vector
%            x    = explanatory variables matrix
%          Symat  = matrix from mess_g3 
%                   containing a mesh of Sy for alph,rho,m values
%            beta = kx1 current bhat vector
%            rho  = current rho value (for lookup)
%           neigh = current m value (for lookup)
%            amat = matrix of alpha values in the mesh
%            rmat = matrix of rho values in the mesh
%            mmat = matrix of neigh values in the mesh
% ---------------------------------------------------
%  RETURNS: a conditional used in Metropolis-Hastings sampling
%  NOTE: called only by mess_g3
%  --------------------------------------------------
%  SEE ALSO: mess_g3, mess_g3d
% ---------------------------------------------------

% written by: James P. LeSage 1/2000
% University of Toledo
% Department of Economics
% Toledo, OH 43606
% [e-mail]


n = length(ys);
% lookup Sy based on alpha, rho values
gsize = rmat(2,1) - rmat(1,1);
i1 = find(rmat <= rho + gsize);
i2 = find(rmat <= rho - gsize);
i1 = max(i1);
i2 = max(i2);
indexr = round((i1+i2)/2);
if isempty(indexr)
indexr = 1;
end;
gsize = mmat(2,1) - mmat(1,1);
i1 = find(mmat <= neigh + gsize);
i2 = find(mmat <= neigh - gsize);
i1 = max(i1);
i2 = max(i2);
indexm = round((i1+i2)/2);
if isempty(indexm)
indexm = 1;
end;

Ycap = squeeze(Ymat(:,:,indexr,indexm));

[junk nq] = size(Ycap);
nq1 = nq-1;
v = ones(nq,1);
for i=2:nq;
v(i,1) = alpha.^(i-1);
end;
W = (1./[1 cumprod(1:nq1)]);
Sy = Ycap*diag(W)*v;
e = Sy - xs*beta;
epe = e'*e;
cout = -(n/2)*log(epe); 


function cout = cn_mess3(alpha,ys,xs,Ymat,beta,sige,rho,neigh,rmat,mmat,a,B);
% PURPOSE: evaluate the conditional distribution of alpha 
%          for the Bayesian mess_g3 model
% ---------------------------------------------------
%  USAGE: cout = c_mess3(alpha,y,x,Symat,beta,sige,rho,neigh,rmat,mmat,a,B)
%  where:  alpha  = matrix exponential alpha spatial parameter
%            y    = dependent variable vector
%            x    = explanatory variables matrix
%          Symat  = matrix from mess_g3 
%                   containing a mesh of Sy for rho,m values
%            beta = kx1 current bhat vector
%            sige = 1x1 current sige value
%            rho  = current rho value (for lookup)
%           neigh = current m value (for lookup)
%            rmat = matrix of rho values in the mesh
%            mmat = matrix of neigh values in the mesh
%               a = prior mean for alpha     (optional input)
%               B = prior variance for alpha (optional input)
% ---------------------------------------------------
%  RETURNS: a conditional used in Metropolis-Hastings sampling
%  NOTE: called only by mess_g3
%  --------------------------------------------------
%  SEE ALSO: mess_g3, mess_g3d
% ---------------------------------------------------

% written by: James P. LeSage 1/2000
% University of Toledo
% Department of Economics
% Toledo, OH 43606
% [e-mail]


n = length(ys);
% lookup Sy based on alpha, rho values
gsize = rmat(2,1) - rmat(1,1);
i1 = find(rmat <= rho + gsize);
i2 = find(rmat <= rho - gsize);
i1 = max(i1);
i2 = max(i2);
indexr = round((i1+i2)/2);
if isempty(indexr)
indexr = 1;
end;
gsize = mmat(2,1) - mmat(1,1);
i1 = find(mmat <= neigh + gsize);
i2 = find(mmat <= neigh - gsize);
i1 = max(i1);
i2 = max(i2);
indexm = round((i1+i2)/2);
if isempty(indexm)
indexm = 1;
end;

Ycap = squeeze(Ymat(:,:,indexr,indexm));

[junk nq] = size(Ycap);
nq1 = nq-1;
v = ones(nq,1);
for i=2:nq;
v(i,1) = alpha.^(i-1);
end;
W = (1./[1 cumprod(1:nq1)]);
Sy = Ycap*diag(W)*v;
e = Sy - xs*beta;
epe = e'*e;

if nargin == 10
cout =  -epe/(2*sige);
elseif nargin == 12
B = B*sige;
cout = -epe/(2*sige) - 0.5*(((alpha-a)^2)/B);
else
error('cn_mess3: Wrong # of inputs arguments');
end;



function cout = nn_mess3(neigh,ys,xs,Ymat,beta,alpha,rho,rmat,mmat);
% PURPOSE: evaluate the conditional distribution of alpha 
%          for the Bayesian mess_g2 model
% ---------------------------------------------------
%  USAGE: cout = c_mess3(alpha,y,x,Symat,beta,alpha,rho,neigh,amat,rmat,mmat)
%  where:  alpha  = matrix exponential alpha spatial parameter
%            y    = dependent variable vector
%            x    = explanatory variables matrix
%          Symat  = matrix from mess_g3 
%                   containing a mesh of Sy for rho,m values
%            beta = kx1 current bhat vector
%            rho  = current rho value (for lookup)
%           neigh = current m value (for lookup)
%            rmat = matrix of rho values in the mesh
%            mmat = matrix of neigh values in the mesh
% ---------------------------------------------------
%  RETURNS: a conditional used in Metropolis-Hastings sampling
%  NOTE: called only by mess_g3
%  --------------------------------------------------
%  SEE ALSO: mess_g3, mess_g3d
% ---------------------------------------------------

% written by: James P. LeSage 1/2000
% University of Toledo
% Department of Economics
% Toledo, OH 43606
% [e-mail]


n = length(ys);
% lookup Sy based on alpha, rho values
gsize = rmat(2,1) - rmat(1,1);
i1 = find(rmat <= rho + gsize);
i2 = find(rmat <= rho - gsize);
i1 = max(i1);
i2 = max(i2);
indexr = round((i1+i2)/2);
if isempty(indexr)
indexr = 1;
end;
gsize = mmat(2,1) - mmat(1,1);
i1 = find(mmat <= neigh + gsize);
i2 = find(mmat <= neigh - gsize);
i1 = max(i1);
i2 = max(i2);
indexm = round((i1+i2)/2);
if isempty(indexm)
indexm = 1;
end;

Ycap = squeeze(Ymat(:,:,indexr,indexm));

[junk nq] = size(Ycap);
nq1 = nq-1;
v = ones(nq,1);
for i=2:nq;
v(i,1) = alpha.^(i-1);
end;
W = (1./[1 cumprod(1:nq1)]);
Sy = Ycap*diag(W)*v;
e = Sy - xs*beta;
epe = e'*e;
cout = -(n/2)*log(epe);