*
* Example 8.8.7 from pp 306-309
*
open data goals.tsm
data(format=free,org=columns) 1 57 england year
*
* This extracts the subsample where there is valid (non-missing) data for England.
*
sample(smpl=%valid(england)) england / scores
stats scores
*
* A "static" Poisson model can be estimated by using DDV (Discrete Dependent
* Variables) with TYPE=COUNT. With a fixed probability, the only explanatory
* variable is the constant. Note that this shows the overdispersion by the
* standardized variance, which should be 1 if a Poisson is actually correct.
*
ddv(type=count) scores
# constant
*
* alpha and lambda will be generated recursively during evaluations of the
* likelihood function.
*
set alpha  = scores
set lambda = 0.0
*
nonlin delta
*
* We have to give an initial guess value for delta because at the default (zero),
* the likelihood function isn't defined.
*
compute delta=.500
*
* This is the main formula for updating the two "state" variables, and evaluating
* the log likelihood. Note that there's a typo in the book: the final parameter
* in the nb density should be lambda/(1+lambda).
*
frml logl = lambda=delta+delta*lambda{1},alpha=delta*scores+delta*alpha{1},$
    log(%negbin(scores,alpha{1},lambda{1}/(1+lambda{1})))
*
* This sets the pre-sample values for lambda and delta. A "startup" formula like
* this is executed just once per function evaluation. This is the sequence which
* generates the estimates in the data. You generated this by starting at lambda
* and alpha=0, then updating for the first two data points (with 0 and 1 goal
* respectively), with the estimation starting at period 3.
*
frml init = lambda{1}=delta*(1+delta),alpha{1}=delta*scores{1}
*
* The likelihood function is maximized over the range 3-52. You can't start
* earlier than 3, because 2 is the first period with a non-zero value for "score".
*
maximize(startup=init,method=bfgs) logl 3 52
*
* Compute the expected scores over the sample period
*
set escore = alpha/lambda
graph(footer="Figure 8-10 One-step predictors of goal data",overlay=dots,ovsame) 2
# escore
# scores
*
* This will compute the predicted probabilities given the parameters at the end
* of the period.
*
dec vect predicted(6)
ewise predicted(i)=%negbin(i-1,alpha(52),lambda(52)/(1+lambda(52)))
disp predicted