* * ARMAGIBBS.RPF * Example of analysis of an ARMA model using Gibbs sampling. * This uses a "flat" prior on the coefficients other than restricting * them to the interval (-1,1). * open data agquarterly.xls calendar(q) 1960:1 data(format=xls,org=columns) 1960:1 2002:1 ppi graph(header="Producer Price Index, 1995=100") # ppi * set logppi = log(ppi) * boxjenk(diffs=1,constant,ar=1,ma=1,maxl) logppi * * The proposal density will be the asymptotic distribution from maximum * likelihood, thus, the mean is %BETA and the factor of the covariance * matrix is %SEESQ*%XX. We're also starting the sampling from these * estimates, so LOGPLAST is %LOGL (we're using a non-informative prior) * and LOGQLAST is 0. * compute bboxj=%beta compute fxx =%decomp(%seesq*%xx) compute logplast=%logl compute logqlast=0.0 compute bdraw=bboxj * compute nburn =2500 compute ndraws=10000 * dec series[vect] bgibbs gset bgibbs 1 ndraws = %zeros(%nreg,1) * compute accept=0 infobox(action=define,lower=-nburn,upper=ndraws,progress) "Independence Metropolis" do draw=-nburn,ndraws * * Draw the proposal and save the log kernel for the draw * compute [vect] btest=bboxj+%ranmvt(fxx,10.0) compute logqtest=%ranlogkernel() * * Reject if AR or MA coefficient is out of range. * if abs(btest(2))>=1.0.or.abs(btest(3))>=1.0 goto reject * * Evaluate the likelihood and do the Metropolis test. * boxjenk(diffs=1,constant,ar=1,ma=1,maxl,initial=btest,method=evaluate) logppi compute logptest=%logl compute alpha =exp(logptest-logplast-logqtest+logqlast) if alpha>1.0.or.%uniform(0.0,1.0)