*
* @KSCPostDraw( options ) y zterms zsqterms zdraw
*
* KSCPostDraw uses the rejection method to draw from the posterior
* formed by multiplying a Normal by a log inverse gamma. This comes up
* in the analysis of the stochastic volatility model. The data evidence
* on the volatility for a given data point takes the form of the log
* inverse gamma, that is, the log volatility is an inverse gamma. The
* information from neighboring data points is Gaussian.
*
* Parameters:
*   y         observed data
*
*   zterms    the linear and quadratic terms in the log density for
*   zsqterms   the Gaussian distribution from the neighbors
*
*   zdraw     (output) draw from posterior
*
* Options:
*   ITERS=number of Newton's method steps to refine expansion point [5]
*     0 is the original proposal of Kim, Shephard and Chib. That can fail
*     to produce an adequate approximation. ITERS=5 seems to work well
*     in almost all cases.
*
* References:
*  Kim, Shephard and Chib(1998), "Stochastic Volatility: Likelihood
*    Inference and Comparison with ARCH Models", Review of Economic
*    Studies, vol 65, pp 361-93
*
procedure KSCPostDraw y zterms zsqterms zdraw
type real y zterms zsqterms *zdraw
*
option integer iters 5
option switch  trace 0
*
local real z0 zadjust zsqadjust zvar zmean zvar_s zmean_s loga
local integer it trys
dec integer failed
*
* Mean and variance of the Gaussian component
*
compute zmean_s=zterms/zsqterms,zvar_s=1.0/zsqterms
*
* Kim-Shephard-Chib rejection sampling with initial quadratic expansion
* to approximate the posterior mode. ITERS=0 gives original KSC method.
*
compute z0=zmean_s
do it=1,iters
    compute zsqadjust=y^2*exp(-z0)*.5
    compute zadjust  =-.50+.50*y^2*exp(-z0)*(1+z0)
    compute zvar     =1.0/(zsqterms+zsqadjust)
    compute zmean    =(zterms+zadjust)*zvar
    if abs(z0-zmean)/sqrt(zvar)<.001
       break
    compute z0=zmean
end do iters
*
*  KSC rejection - linearize around expansion point
*
compute zsqadjust=0.0
compute zadjust  =-.50+.50*y^2*exp(-z0)
compute zvar     =1.0/(zsqterms+zsqadjust)
compute zmean    =(zterms+zadjust)*zvar
do trys=1,1000
    compute zdraw=zmean+%ran(sqrt(zvar))
    compute loga=-.50*y^2*(exp(-zdraw)-exp(-z0)+exp(-z0)*(zdraw-z0))
    if trace
       disp trys zmean_s zmean zvar_s zvar zdraw exp(loga)
    if %uniform(0.0,1.0)<exp(loga)
       break
end do trys
compute failed=0
if trys==1000 {
   disp "Failed on" y zterms zsqterms
   compute failed=1
}
*compute totaltrys=totaltrys+trys
end
