/* ===========================================================
FIND MAXIMUM LIKELIHOOD ESTIMATES, IF DESIRED
=============================================================*/
kc = 1;
if kmle == 1; @ this starts the MLE section @
/* ------------------------------------------------------------------------------------------------------------------------
=== SET STARTING VALUES FOR NUMERICAL OPTIMIZATION === */
th0 = ( 1./(sigx*sqrt(k)) ) | 0.5;
th0';
@ thx = 0.084 | 0.001 | 1.85 ; h = eye(k+1); i = 1; psi = zeros(1,k+2); goto jump4; @
proc startval; @ This defines starting value for iteration to be th @
retp(th0); endp;
/* -----------------------------------------------------------------------------------------------------------------------------
=== SET GLOBAL CONTROL PARAMETERS FOR NUMERICAL
OPTIMIZATION === */
#include optmum.ext;
_opalgr = 2; @ This chooses BFGS optimization @
_opmiter = 150; @ This controls the maximum number of iterations @
_opshess = (1.e-4)*eye(k+1);
/* -----------------------------------------------------------------------------------------------------------------------------
=== CALL GAUSS PROCEDURES FOR NUMERICAL OPTIMIZATION === */
kc=1;
{thx,f,g,h} =optmum(&conc,startval);
"";"";"Final results";"MLE as parameterized for numerical optimization ";
"Coefficients:";thx';
"";"Value of log likelihood:";;-f;
"";"Gradient vector:";g';
"==================================";
kc = 2;
call conc(thx);
/* ------------------------------------------------------------------------------------------------------------------------
=== CALCULATE STANDARD ERRORS === */
kc = 3;
ks = 0;
thmore = conc(thx);
kc = 1;
g = gradp(&full,thx|thmore);
"";"Numerically calculated gradient:";g;
"";"";
h = (hessp(&full,thx|thmore));
va = eigrs(h);
if minc(eigrs(h)) <= 0;
"Negative of Hessian is not positive definite";
"Either you have not found local maximum, or else estimates are up "
"against boundary condition. In latter case, impose the restricted "
"params rather than estimate them to calculate standard errors";
else;
h = invpd(h);
std = diag(h)^.5;
"For vector of coefficients parameterized as follows,";(thx | thmore)';
"the standard errors are";std';
endif;
"variance covariance matrix is";
format /m3;
h;
format /m1;
h;
/* -------------------------------------------------------------------------------------------------------------------------
=== GRAPH FUNCTIONS OF INTEREST === */
_pgrid = 3 | 0;
if kex == 1;
@ plot figure 4 for MLE @
_ptek = "graph4ml.tkf";
title("Figure 4 (MLE)");
call graph1dim(1,2,xbar);
chis = conc(thx);
y1hat = 0.6*xs[.,1].*(xs[.,1].>zeros(ks,1)) + 0.2*xbar[2,1]*ones(ks,1);
xy(xs[.,1],chis~y1hat);
@ plot figure 5 for MLE @
_ptek = "graph5ml.tkf";
title("Figure 5 (MLE)");
call graph1dim(2,2,xbar);
chis = conc(thx);
xy(xs[.,2],chis);
elseif kex == 2;
/* _pgrid = 0 | 0;
_ptek = "graph6.tkf";
title("Figure 6");
{x3,y3,z3} = graph2dim(1,2,2,thx,xbar);
contour(x3',y3,z3); */
elseif kex == 3;
_ptek = "graph7ml.tkf";
title("Figure 7 (MLE)");
call graph1dim(3,2,xbar);
chis = conc(thx);
xy(xs[.,3],chis);
endif;
if kex == 3;
@ This section calculates values used in Figure 8 @
kc = 1;
ks = 61;
"Conditional mean at different dates";
xs = seqa(xbar[1,1]-2*sigx[1,1],(4*sigx[1,1])/(ks-1),ks)~
(ones(ks,1)*(-0.7))~(ones(ks,1)*1955);
chis = conc(thx);
xs~chis;
"";
xs = seqa(xbar[1,1]-2*sigx[1,1],(4*sigx[1,1])/(ks-1),ks)~
(ones(ks,1)*(12.3))~(ones(ks,1)*1975);
chis = conc(thx);
xs~chis;
"";
xs = seqa(xbar[1,1]-2*sigx[1,1],(4*sigx[1,1])/(ks-1),ks)~
(ones(ks,1)*(3.9)~ones(ks,1)*1985);
chis = conc(thx);
xs~chis;
endif;
endif; @ this ends the MLE section @
/* ============================================================
PERFORM IMPORTANCE SAMPLING, IF DESIRED
=============================================================*/
/* ============================================================
SET INITIAL VALUES FOR PARAMETERS
============================================================ */
jump4:
@ ---------------------- parameters of prior for theta ------------------------------------------ @
meta = 0; @ meta is mean of log of eta @
seta = 1.0; @ seta is std deviation of log of eta @
mg = ln(1./(sigx*sqrt(k))); @ mg is k x 1 vector of prior means for ln(g) @
taug = 1.0*ones(k,1); @ taug is prior std. deviation for ln(g) @
mth = mg | meta;
sigth = taug | seta;
@ ---------------------- parameters of prior for psi ------------------------------------------------ @
nsig = 0.25; @ prior is (1/sigma^2) ~ gamma(nsig, lamsig) @
lamsig = nsig*0.5*(sigy^2);
hbet = n;
mbet = zeros(klin,1);
mbet[1,1] = meanc(y);
@ prior is beta ~ N(mbet,hbet*sigma^2*invpd(xwhole'*xwhole)) @
@ ----------parameters that determine importance sampling density -------------- @
montdf = 2; @ montdf is degrees of freedom for student t
generated importance density @
montfact = sqrt(2); @ montfact is factor by which standard deviations are
increased to obtain std. deviation of importance density @
pprob = 0.5; @ pprob = probability of drawing from the Student t versus
the spread-out prior @
"";"prior distributions:";
" ln(eta) ~ Normal(";;meta;;",";;seta;"^2)";
ii = 1;
do until ii > k;
" ln(g(";;ii;;")) ~ Normal(";;mg[ii,1];;",";;taug[ii,1];;"^2)";
ii = ii+1;
endo;
" (1/sigma^2) | eta ~ Gamma(";;nsig;;",";; lamsig;;")";
" beta | sigma,theta ~ N(mbet,hbet*sigma*invpd(xwhole'*xwhole)) ";
" mbet:";;mbet';
" hbet:";;hbet;
"";
"Importance density:";
" degreees of freedom of t distribution";;montdf;
" fraction of observations from t distribution";;pprob;
" factor by which std. deviation of importance density exceeds prior";;montfact;
@ ----------------------- parameters to control monte carlo runs -----------------------@
nmonte = 20000; @ nmonte is number of monte carlo draws generated @
"Number of monte carlo draws";;nmonte;
"";
if kmle == 1;
thmle = abs(thx);
pmix = h[1:k+1,1:k+1];
else;
"you must input a mean and variance for the student t mixing distribution";
endif;
kmle = 2;
pmix = pmix*montfact^2;
uxmix = zeros(nmonte,1); @uxmix can be used to keep track of
which component of importance density was used @
@ -------- calculation of constant terms for importance density ------------------- @
tprec = invpd(pmix);
tdet = detl;
tc = gamma( (montdf+k+1)/2 )/gamma(montdf/2);
tc = tc / ( sqrt(tdet) * (montdf * 3.14159)^( (k+1) / 2) );
nc =1./( sqrt(2*3.14159).*sigth); @ note that nc is (k +1) x 1 vector @
lnnc =ln(nc);
pmix = chol(pmix);
/* =============================================================
INCLUDE NECESSARY PROCEDURES
============================================================= */
psi = 0; i = 0; @ needed to avoid tripping GAUSS error compiling code @
#include bayproc2;
proc normgam(sigtry,bettry,ngam,lamgam,betsy,varsy);
@ this proc evaluates the log of the product of a gamma(ngam,lamgam) density
for (1/sigtry) and a N(betsy,sigtry*varsy) density for bettry at the
points sigtry and bettry @
local f1,f2,m1,m2;
if ngam > 10;
f1 = ngam*ln(lamgam) - (ngam - 1)*ln(sigtry) - (lamgam/sigtry) - lnfact(ngam-1);
else;
f1 = ngam*ln(lamgam) - (ngam - 1)*ln(sigtry) - (lamgam/sigtry) - ln(gamma(ngam));
endif;
m1 = invpd(varsy);
m2 = detl;
f2 = (-klin/2)*ln(2*3.1415927*sigtry) -(1/2)*ln(m2);
f2 = f2 - (1/(2*sigtry))*(bettry-betsy)'*m1*(bettry-betsy);
retp(f1 + f2);
endp;
proc priadj(th);
@ this proc adds log of prior p(theta) to log of f(y|X,theta) to
arrive at log of f(y,theta|X) @
local q;
if ndpchk(16);
"underflow before calling priadj";
endif;
q = -((ln(th) - mth)^2)./(2*(sigth^2));
q = q + lnnc - ln(th);
q = sumc(q);
if kc > 2;
"log of prior is";;q;
endif;
q = q + baby(th);
if ndpchk(16);
"underflow after calling priadj";
endif;
retp(q);
endp;
发帖
雷达卡
京公网安备 11010802022788号
