*
* Illustration from section 7.7.8, pp 330-334
*
open data d-ibmln98wm.txt
data(format=free,org=columns) 1 9190 time ibmcret
set ibmcret = ibmcret*100.0
*
open data d-ibml25x.txt
data(format=free,org=columns) 1 9190 x1 x2 x3 x4 x5
*
set flip = -ibmcret
*
nonlin(parmset=betaparms)  b0 b1 b2 b3 b4 b5
nonlin(parmset=kparms)     k0 k1 k2 k3 k4 k5
nonlin(parmset=alphaparms) a0 a1 a2 a3 a4 a5
*
compute eta=2.50
frml bf = b0+b1*x1+b2*x2+b3*x3+b4*x4+b5*x5
frml kf = k0+k1*x1+k2*x2+k3*x3+k4*x4+k5*x5
frml af = a0+a1*x1+a2*x2+a3*x3+a4*x4+a5*x5
*
frml lambdaf = (9190.0/252.0)*%max((1-kf*(eta-bf)/exp(af)),0)**(1.0/kf)
*
compute k0=-.143
compute b0=2.00
compute a0=0.00
*
compute b1=b2=b3=b4=b5=0.0
compute k1=k2=k3=k4=k5=0.0
compute a1=a2=a3=a4=a5=0.0
*
find(method=bfgs,parmset=betaparms+alphaparms+kparms,stderrs) maximum %logl
   sstats / lambdaf>>ssum
   sstats(smpl=(flip>eta)) / %loggpdensity(flip,kf,eta,exp(af))+log(lambdaf)>>esum
   compute %logl=-ssum/252.0+esum
end
*
nonlin(parmset=drops) b1=0 b2=0 b3=0 a1=0 a2=0 k1=0 k2=0
find(method=bfgs,parmset=betaparms+alphaparms+kparms+drops,stderrs) maximum %logl
   sstats / lambdaf>>ssum
   sstats(smpl=(flip>eta)) / %loggpdensity(flip,kf,eta,exp(af))+log(lambdaf)>>esum
   compute %logl=-ssum/252.0+esum
end