*
* Example 7.7 from pp 323-329
*
open data d-ibmln98.dat
data(format=free,org=columns) 1 9190 ibmlog
set flip   = -ibmlog
*
spgraph(vfields=2,footer="Figure 7.6 Plots for daily negative IBM log returns")
set smpl = flip>=0.025
@qplot(smpl=smpl,distrib=expon,bound=.025,$
  title="(a) Q Plot of negative daily returns vs exponential") flip
@meplot(title="(b) Mean excess plot") flip
spgraph(done)
*
nonlin k alpha beta
compute k=-.143
compute beta=2.00
compute alpha=1.00
*
dec real eta n0
*
frml bipoisson = lambda=(9190.0/252.0)*%max((1-k*(eta-beta)/alpha),0)**(1.0/k),$
        %loggpdensity(flip,k,eta,alpha)+log(lambda)-lambda/n0
dofor [real] eta = 3.00 2.50 2.00
   sstats / (flip>eta)>>n0
   maximize(smpl=flip>eta) bipoisson
   disp "VaR(5%)" 100000*%invgev((1-.05)**252,k,beta,alpha)  $
        "VaR(1%)" 100000*%invgev((1-.01)**252,k,beta,alpha)
end dofor
*
* Same analysis with means extracted
*
diff(center) flip / cflip
frml bipoisson = lambda=(9190.0/252.0)*%max((1-k*(eta-beta)/alpha),0)**(1.0/k),$
         %loggpdensity(cflip,k,eta,alpha)+log(lambda)-lambda/n0
*
dofor [real] eta = 3.00 2.50 2.00
   sstats / (cflip>eta)>>n0
   maximize(smpl=cflip>eta) bipoisson
   disp "VaR(5%)" 100000*%invgev((1-.05)**252,k,beta,alpha)  $
        "VaR(1%)" 100000*%invgev((1-.01)**252,k,beta,alpha)
end dofor
*
*  Switch back to the model in terms of uncentered returns
*
compute eta=2.50
frml bipoisson = lambda=(9190.0/252.0)*%max((1-k*(eta-beta)/alpha),0)**(1.0/k),$
        %loggpdensity(flip,k,eta,alpha)+log(lambda)-lambda/n0
sstats / (flip>eta)>>n0
maximize(smpl=flip>eta) bipoisson
compute var=%invgp(1-.01*9190/n0,k,eta,alpha)
compute es =(var+alpha+k*eta)/(1+k)
disp "99%-ile"	"VAR=" var "ES=" es