prt_gibbs.m

时间序列分析中常用到的matlab代码

end;
end;

fprintf(fid,'\n');
fprintf(fid,'Bayesian Model Averaging Estimates \n');
if nflag == 1
fprintf(fid,'Dependent Variable = %16s \n',vnames(1,:));
end;
fprintf(fid,'R-squared      = %9.4f \n',results.rsqr);
fprintf(fid,'sigma^2        = %9.4f \n',results.sige);
fprintf(fid,'Nobs, Nvars    = %6d,%6d \n',results.nobs,results.nvar);
fprintf(fid,'ndraws         = %6d \n',results.ndraw);
fprintf(fid,'nu,lam,phi     = %6.3f,%6.3f,%6.3f \n',results.nu,results.lam,results.phi);
fprintf(fid,'# of models    = %6d   \n',results.nmod);
fprintf(fid,'time(seconds)  = %9.4f \n',results.time);

fprintf(fid,'***************************************************************\n');

outi = find(results.prob > 0.01); % limit printing to models > 1 percent post prob
nmod = length(outi);
out = [results.model(outi,:) results.prob(outi,1)*100 results.visit(outi,1)];
 fmt = '%5d';
 cnames = Vname{1};
 for i=2:results.nvar
 cnames = strvcat(cnames,Vname{i});
 fmt = strvcat(fmt,'%5d');
 end;
 fmt = strvcat(fmt,'%6.3f');   
 fmt = strvcat(fmt,'%5d');   
cnames = strvcat(cnames,'Prob','Visit');

rnames = 'Model';
for i=1:nmod;
    rnames = strvcat(rnames,['model ' num2str(i)]); 
end;

min.cnames = cnames;
min.rnames = rnames;
min.fid = fid;
min.fmt = fmt;
fprintf(fid,'Model averaging information \n');
mprint(out,min);
fprintf(fid,'***************************************************************\n');
fprintf(fid,'      Posterior Estimates \n');
tstat = results.tstat;
bhat = results.beta;
nvar = results.nvar+1; % need to include the constant term

% now print coefficient estimates, t-statistics and probabilities

% find t-stat marginal probabilities
tout = tdis_prb(tstat,results.nobs);

% column labels for printing results
vstring = 'Variable';
bstring = 'Coefficient';
tstring = 't-statistic';
pstring = 't-probability';

tmp = [bhat tstat tout];

cnames = strvcat(bstring,tstring,pstring);
rnames = vstring;
if nflag == 0
 rnames = strvcat(rnames,'const');
 for i=1:results.nvar
 rnames = strvcat(rnames,Vname{i});
 end;
else
 rnames = strvcat(rnames,'const');
 for i=1:results.nvar
  rnames = strvcat(rnames,Vname{i});
 end;
end;

in2.fmt = '%16.6f';
in2.fid = fid;
in2.rnames = rnames;
in2.cnames = cnames;

mprint(tmp,in2);

case {'probit_g'} % <=================== heteroscedastic probit model    
    
% we handle these differently depending on the model
if ( nflag == 0) %  no variable names supplied, make some up
Vname = [];
for i=1:nvar
    Vname{i} = str2mat(['variable   ',num2str(i)]);
end;
end;

bhat = mean(results.bdraw);  % calculate means and std deviations
bhat = bhat';
bstd = std(results.bdraw);
bstd = bstd';

if strcmp(results.pflag,'tstat')
 tstat = bhat./bstd;
 % find t-stat marginal probabilities
 tout = tdis_prb(tstat,results.nobs);
else % find plevels
 for i=1:results.nvar;
 if bhat(i,1) > 0
 cnt = find(results.bdraw(:,i) > 0);
 tout(i,1) = 1 - (length(cnt)/(results.ndraw-results.nomit));
 else
 cnt = find(results.bdraw(:,i) < 0);
 tout(i,1) = 1 - (length(cnt)/(results.ndraw-results.nomit));
 end; % end of if - else
 end; % end of for loop
end; 

nobs = results.nobs;
nvar = results.nvar;
y = results.y;
zip = length(find(y == 0));
one = nobs - zip;
fprintf(fid,'\n');
fprintf(fid,'Bayesian Heteroscedastic Probit Model Gibbs Estimates \n');
if nflag == 1
fprintf(fid,'Dependent Variable = %16s \n',vnames(1,:));
end;
fprintf(fid,'McFadden R^2    = %9.4f \n',results.r2mf);
fprintf(fid,'Estrella R^2    = %9.4f \n',results.rsqr);
fprintf(fid,'Nobs, Nvars     = %6d,%6d \n',results.nobs,results.nvar);
fprintf(fid,'# 0, 1 y-values = %6d,%6d \n',zip,nobs-zip);
fprintf(fid,'ndraws,nomit    = %6d,%6d \n',results.ndraw,results.nomit);
fprintf(fid,'time in secs    = %9.4f\n',results.time);
rmean = mean(results.rdraw);
if rmean ~= 0
fprintf(fid,'rmean           = %9.4f \n',rmean);
else
fprintf(fid,'r-value         = %6d  \n',results.r);
end;

fprintf(fid,'***************************************************************\n');

vstring = 'Variable';
bstring = 'Prior Mean';
tstring = 'Std Deviation';

tmp = [results.pmean results.pstd];

cnames = strvcat(bstring,tstring);
rnames = vstring;
for i=1:nvar
rnames = strvcat(rnames,Vname{i});
end;
pin.fmt = '%16.6f';
pin.fid = fid;
pin.cnames = cnames;
pin.rnames = rnames;

mprint(tmp,pin);
fprintf(fid,'***************************************************************\n');
fprintf(fid,'      Posterior Estimates \n');
% now print coefficient estimates, t-statistics and probabilities


% column labels for printing results
vstring = 'Variable';
bstring = 'Coefficient';

if strcmp(results.pflag,'tstat') % depends on pflag argument
tstring = 't-statistic';
pstring = 't-probability';
tmp = [bhat tstat tout];
else
tstring = 'Std Deviation';
pstring = 'p-level';
tmp = [bhat bstd tout];
end;


cnames = strvcat(bstring,tstring,pstring);
rnames = vstring;
for i=1:nvar
rnames = strvcat(rnames,Vname{i});
end;
in.fmt = '%16.6f';
in.fid = fid;
in.cnames = cnames;
in.rnames = rnames;
mprint(tmp,in);

case {'tobit_g'} % <=================== heteroscedastic tobit model    
    
% we handle these differently depending on the model
if ( nflag == 0) %  no variable names supplied, make some up
Vname = [];
for i=1:nvar
    Vname{i} = str2mat(['variable   ',num2str(i)]);
end;
end;

yact = results.y;
bhat = mean(results.bdraw);  % calculate means and std deviations
bhat = bhat';
bstd = std(results.bdraw);
bstd = bstd';
if strcmp(results.pflag,'tstat')
 tstat = bhat./bstd;
 % find t-stat marginal probabilities
 tout = tdis_prb(tstat,results.nobs);
else % find plevels
 for i=1:results.nvar;
 if bhat(i,1) > 0
 cnt = find(results.bdraw(:,i) > 0);
 tout(i,1) = 1 - (length(cnt)/(results.ndraw-results.nomit));
 else
 cnt = find(results.bdraw(:,i) < 0);
 tout(i,1) = 1 - (length(cnt)/(results.ndraw-results.nomit));
 end; % end of if - else
 end; % end of for loop
end; 

nobs = results.nobs;
nvar = results.nvar;
yhat = results.x*bhat;
resid = yact - yhat;
sigu = resid'*resid;
ym = yact - ones(nobs,1)*mean(yact);
rsqr1 = sigu/(nobs-nvar);
rsqr2 = ym'*ym;
rsqr = 1.0 - rsqr1/rsqr2; % conventional r-squared

fprintf(fid,'\n');
fprintf(fid,'Bayesian Heteroscedastic Tobit Model Gibbs Estimates \n');
if nflag == 1
fprintf(fid,'Dependent Variable = %16s \n',vnames(1,:));
end;
fprintf(fid,'R-squared      = %9.4f \n',rsqr); % based on draws
fprintf(fid,'sigma^2        = %9.4f \n',mean(results.sdraw));
fprintf(fid,'nu,d0          = %6d,%6d \n',results.nu,results.d0);
fprintf(fid,'Nobs, Nvars    = %6d,%6d \n',results.nobs,results.nvar);
fprintf(fid,'# of censored  = %6d \n',results.nobsc);
fprintf(fid,'ndraws,nomit   = %6d,%6d \n',results.ndraw,results.nomit);
fprintf(fid,'time in secs   = %9.4f\n',results.time);
rmean = mean(results.rdraw);
if rmean ~= 0
fprintf(fid,'rmean          = %9.4f \n',rmean);
else
fprintf(fid,'r-value        = %6d  \n',results.r);
end;

fprintf(fid,'***************************************************************\n');

vstring = 'Variable';
bstring = 'Prior Mean';
tstring = 'Std Deviation';

tmp = [results.pmean results.pstd];

cnames = strvcat(bstring,tstring);
rnames = vstring;
for i=1:nvar
rnames = strvcat(rnames,Vname{i});
end;
pin.fmt = '%16.6f';
pin.fid = fid;
pin.cnames = cnames;
pin.rnames = rnames;

mprint(tmp,pin);
fprintf(fid,'***************************************************************\n');
fprintf(fid,'      Posterior Estimates \n');
% now print coefficient estimates, t-statistics and probabilities

% column labels for printing results
vstring = 'Variable';
bstring = 'Coefficient';
if strcmp(results.pflag,'tstat') % depends on pflag argument
tstring = 't-statistic';
pstring = 't-probability';
tmp = [bhat tstat tout];
else
tstring = 'Std Deviation';
pstring = 'p-level';
tmp = [bhat bstd tout];
end;


cnames = strvcat(bstring,tstring,pstring);
rnames = vstring;
for i=1:nvar
rnames = strvcat(rnames,Vname{i});
end;
in.fmt = '%16.6f';
in.fid = fid;
in.cnames = cnames;
in.rnames = rnames;
mprint(tmp,in);


otherwise
error('results structure not known by prt_gibbs function');

end;