semp_g.m

计量工具箱

%          ldetflag,eflag,order,iter,novi_flag,c,T,prior_beta,cc,metflag,a1,a2,inform_flag] = 
%                           sem_parse(prior,k)
% where info contains the structure variable with inputs 
% and the outputs are either user-inputs or default values
% ---------------------------------------------------

% set defaults

eflag = 1;     % default to not computing eigenvalues
ldetflag = 1;  % default to 1999 Pace and Barry MC determinant approx
order = 50;    % there are parameters used by the MC det approx
iter = 30;     % defaults based on Pace and Barry recommendation
rmin = -1;     % use -1,1 rho interval as default
rmax = 1;
detval = 0;    % just a flag
rho = 0.5;
sige = 1.0;
rval = 4;
mm = 0;
kk = 0;
nu = 0;
d0 = 0;
a1 = 1.01;
a2 = 1.01;
c = zeros(k,1);   % diffuse prior for beta
T = eye(k)*1e+12;
novi_flag = 0;    % default is do vi-estimates
metflag = 1;      % default to Metropolis-Hasting sampling
cc = 0.2;         % initial tuning parameter for M-H sampling
inform_flag = 0;  % flag for diffuse prior on beta


fields = fieldnames(prior);
nf = length(fields);
if nf > 0
 for i=1:nf
    if 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},'eigs')
        eflag = prior.eigs;
    elseif strcmp(fields{i},'dflag')
        metflag = prior.dflag;
    elseif strcmp(fields{i},'a1')
       a1 = prior.a1; 
    elseif strcmp(fields{i},'a2')
       a2 = prior.a2; 
    elseif strcmp(fields{i},'m')
        mm = prior.m;
        kk = prior.k;
        rval = gamm_rnd(1,1,mm,kk);    % initial value for rval   
    elseif strcmp(fields{i},'beta')
        c = prior.beta; inform_flag = 1; % flag for informative prior on beta
    elseif strcmp(fields{i},'bcov')
        T = prior.bcov; inform_flag = 1; % flag for informative prior on beta
    elseif strcmp(fields{i},'rmin')
        rmin = prior.rmin;  
    elseif strcmp(fields{i},'rmax')
        rmax = prior.rmax; 
    elseif strcmp(fields{i},'lndet')
    detval = prior.lndet;
    ldetflag = -1;
    rmin = detval(1,1);
    nr = length(detval);
    rmax = detval(nr,1);
    elseif strcmp(fields{i},'lflag')
        tst = prior.lflag;
        if tst == 0,
        ldetflag = 0; 
        elseif tst == 1,
        ldetflag = 1; 
        elseif tst == 2,
        ldetflag = 2; 
        else
        error('semp_g: unrecognizable lflag value on input');
        end;
    elseif strcmp(fields{i},'order')
        order = prior.order;  
    elseif strcmp(fields{i},'iter')
        iter = prior.iter; 
    elseif strcmp(fields{i},'novi')
        novi_flag = prior.novi;
    end;
 end;
 
else, % the user has input a blank info structure
      % so we use the defaults
end; 


function [rmin,rmax,time2] = sem_eigs(eflag,W,rmin,rmax,n);
% PURPOSE: compute the eigenvalues for the weight matrix
% ---------------------------------------------------
%  USAGE: [rmin,rmax,time2] = far_eigs(eflag,W,rmin,rmax,W)
% where eflag is an input flag, W is the weight matrix
%       rmin,rmax may be used as default outputs
% and the outputs are either user-inputs or default values
% ---------------------------------------------------


if eflag == 0
t0 = clock;
opt.tol = 1e-3; opt.disp = 0;
lambda = eigs(sparse(W),speye(n),1,'SR',opt);  
rmin = 1/lambda;   
rmax = 1;
time2 = etime(clock,t0);
else
time2 = 0;
end;


function [detval,time1] = sem_lndet(ldetflag,W,rmin,rmax,detval,order,iter);
% PURPOSE: compute the log determinant |I_n - rho*W|
% using the user-selected (or default) method
% ---------------------------------------------------
%  USAGE: detval = far_lndet(lflag,W,rmin,rmax)
% where eflag,rmin,rmax,W contains input flags 
% and the outputs are either user-inputs or default values
% ---------------------------------------------------


% do lndet approximation calculations if needed
if ldetflag == 0 % no approximation
t0 = clock;    
out = lndetfull(W,rmin,rmax);
time1 = etime(clock,t0);
tt=rmin:.001:rmax; % interpolate a finer grid
outi = interp1(out.rho,out.lndet,tt','spline');
detval = [tt' outi];
    
elseif ldetflag == 1 % use Pace and Barry, 1999 MC approximation

t0 = clock;    
out = lndetmc(order,iter,W,rmin,rmax);
time1 = etime(clock,t0);
results.limit = [out.rho out.lo95 out.lndet out.up95];
tt=rmin:.001:rmax; % interpolate a finer grid
outi = interp1(out.rho,out.lndet,tt','spline');
detval = [tt' outi];

elseif ldetflag == 2 % use Pace and Barry, 1998 spline interpolation

t0 = clock;
out = lndetint(W,rmin,rmax);
time1 = etime(clock,t0);
tt=rmin:.001:rmax; % interpolate a finer grid
outi = interp1(out.rho,out.lndet,tt','spline');
detval = [tt' outi];

elseif ldetflag == -1 % the user fed down a detval matrix
    time1 = 0;
        % check to see if this is right
        if detval == 0
            error('semp_g: wrgon lndet input argument');
        end;
        [n1,n2] = size(detval);
        if n2 ~= 2
            error('semp_g: wrong sized lndet input argument');
        elseif n1 == 1
            error('semp_g: wrong sized lndet input argument');
        end;          
end;


function  out = sem_marginal(detval,y,x,Wy,Wx,nobs,nvar,a1,a2)
% PURPOSE: returns a vector of the log-marginal over a grid of rho-values
% -------------------------------------------------------------------------
% USAGE: out = sem_marginal(detval,y,x,Wy,Wx,nobs,nvar,a1,a2)
% where:       detval = an ngrid x 2 matrix with rho-values and lndet values
%                  y = y-vector
%                  x = x-matrix
%                 Wy = W*y-vector
%                 Wx = W*x-matrix
%               nobs = # of observations
%               nvar = # of explanatory variables
%                 a1 = parameter for beta prior on rho
%                 a2 = parameter for beta prior on rho
% -------------------------------------------------------------------------
% RETURNS: out = a structure variable
%        out.log = log marginal, a vector the length of detval
%        out.lik = concentrated log-likelihood vector the length of detval
% -------------------------------------------------------------------------
% NOTES: works only for homoscedastic SEM model
% we must feed in ys = sqrt(V)*y, xs = sqrt(V)*X 
% as well as logdetx = log(xs'*xs) for heteroscedastic model
% -------------------------------------------------------------------------

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

nmk = (nobs-nvar)/2;

nrho = length(detval(:,1));
iota = ones(nrho,1);
rvec = detval(:,1);
epe = zeros(nrho,1);
rgrid = detval(1,1)+0.001:0.1:detval(end,1)-0.001;
rgrid = rgrid';

epetmp = zeros(length(rgrid),1);
detxtmp = zeros(length(rgrid),1);
for i=1:length(rgrid);
xs = x - rgrid(i,1)*Wx;
ys = y - rgrid(i,1)*Wy;
bs = (xs'*xs)\(xs'*ys);
e = ys - xs*bs;
epetmp(i,1) = e'*e;
detxtmp(i,1) = det(xs'*xs);
end;

% spline interpolate epetmp
tt=rvec; % interpolate a finer grid
epe = interp1(rgrid,epetmp,rvec,'spline'); 
detx = interp1(rgrid,detxtmp,rvec,'spline'); 


bprior = beta_prior(detval(:,1),a1,a2);
% C is a constant of integration that can vary with nvars, so for model
% comparisions involving different nvars we need to include this
C = log(bprior) + gammaln(nmk) - nmk*log(2*pi) ;
den = detval(:,2) - 0.5*log(detx) - nmk*log(epe);
den = real(den);
out = den + C;

function  out = sem_marginal2(detval,y,x,Wy,Wx,nobs,nvar,a1,a2,c,TI,sige)
% PURPOSE: returns a vector of the log-marginal over a grid of rho-values
%          for the case of an informative prior on beta
% -------------------------------------------------------------------------
% USAGE: out = sem_marginal2(detval,y,x,Wy,Wx,nobs,nvar,a1,a2,c,TI,sige)
% where:       detval = an ngrid x 2 matrix with rho-values and lndet values
%                  y = y-vector
%                  x = x-matrix
%                 Wy = W*y-vector
%                 Wx = W*x-matrix
%               nobs = # of observations
%               nvar = # of explanatory variables
%                 a1 = parameter for beta prior on rho
%                 a2 = parameter for beta prior on rho
% -------------------------------------------------------------------------
% RETURNS: out = a structure variable
%        out.log = log marginal, a vector the length of detval
%        out.lik = concentrated log-likelihood vector the length of detval
% -------------------------------------------------------------------------
% NOTES: works only for homoscedastic SEM model
% we must feed in ys = sqrt(V)*y, xs = sqrt(V)*X 
% as well as logdetx = log(xs'*xs) for heteroscedastic model
% -------------------------------------------------------------------------

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

nmk = (nobs-nvar)/2;

nrho = length(detval(:,1));
iota = ones(nrho,1);
rvec = detval(:,1);
epe = zeros(nrho,1);

rgrid = detval(1,1)+0.001:0.1:detval(end,1)-0.001;
rgrid = rgrid';
epetmp = zeros(length(rgrid),1);
detxtmp = zeros(length(rgrid),1);
Q1 = zeros(length(rgrid),1);
Q2 = zeros(length(rgrid),1);
sTI = sige*TI;

for i=1:length(rgrid);
xs = x - rgrid(i,1)*Wx;
ys = y - rgrid(i,1)*Wy;
bs = (xs'*xs)\(xs'*ys);
beta = (xs'*xs + sTI)\(xs'*ys + sTI*c);
e = ys - xs*bs;
epetmp(i,1) = e'*e;
detxtmp(i,1) = det(xs'*xs);
Q1(i,1) = (c - beta)'*sTI*(c - beta);
Q2(i,1) = (bs - beta)'*(xs'*xs)*(bs - beta);
end;

% spline interpolate epetmp
tt=rvec; % interpolate a finer grid
epe = interp1(rgrid,epetmp,rvec,'spline'); 
detx = interp1(rgrid,detxtmp,rvec,'spline'); 
Q1 = interp1(rgrid,Q1,rvec,'spline'); 
Q2 = interp1(rgrid,Q2,rvec,'spline'); 
bprior = beta_prior(detval(:,1),a1,a2);
% C is a constant of integration that can vary with nvars, so for model
% comparisions involving different nvars we need to include this
C = log(bprior) + gammaln(nmk) - nmk*log(2*pi) ;
den = detval(:,2) - 0.5*log(detx) - nmk*log(epe + Q1 + Q2);
den = real(den);
out = den + C;


function rho = olddraw_rho(detval,y,x,Wy,Wx,V,n,k,rmin,rmax,rho)
% update rho via univariate numerical integration
% for the heteroscedastic model case

nmk = (n-k)/2;
nrho = length(detval(:,1));
iota = ones(nrho,1);
rvec = detval(:,1);
epe = zeros(nrho,1);

for i=1:nrho;
xs = x - rvec(i,1)*Wx;
xs = matmul(xs,sqrt(V));
ys = y - rvec(i,1)*Wy;
ys = ys.*sqrt(V);
bs = (xs'*xs)\(xs'*ys);
e = ys - xs*bs;
epe(i,1) = e'*e;
end;

den = detval(:,2) - nmk*log(epe);
adj = max(den);
den = den - adj;
den = exp(den);


n = length(den);
y = detval(:,1);
x = den;

% trapezoid rule
isum = sum((y(2:n,1) + y(1:n-1,1)).*(x(2:n,1) - x(1:n-1,1))/2);

z = abs(x/isum);
den = cumsum(z);

rnd = unif_rnd(1,0,1)*sum(z);
ind = find(den <= rnd);
idraw = max(ind);
if (idraw > 0 & idraw < nrho)
rho = detval(idraw,1);
end;