/* genoptfit.g  Procedure to obtain fitted call option prices given
                parameters driving the risk neutral jump-diffusion, allowing
                for stochastic volatility and correlation between index and
                unobserved volatility process

           See: David Bates, Review of Financial Studies, v9, n1, 1996 for
                details. The inputs to the moment generating functions were
                translated to Gauss from David Bates Fortan code, which was
                generously supplied by the author.
                      
  FORMAT: cvec=genoptfit(s0,v0,b,tau,strike,theta);


  Inputs:

    S0: scaler or vector of current asset values
    v0: scaler or vector of current instantaneous volatility values
     b: scaler or vector of cost of carry values (see note below)
   tau: scaler or vector time to expiration values in years
strike: scaler or vector of strike prices
 theta: vector of parameters driving the risk neutral process
        "star" parameters represent actual parameters with incorporation
        of risk premiums

         lambdastar=theta[1]  Expected # jumps per year
           kbarstar=theta[2]  Mean % jump, cond'al on jump occuring
              delta=theta[3]  Variance of jump component
              alpha=theta[4]  Longrun Vol. mean
           betastar=theta[5]  Speed of mean reversion
             sigmav=theta[6]  Volatility of Volatility
                rho=theta[7]  Corr. between Underlying and volatility
 __________________________________________________________________
 The actual process is given as follows:

 dS/S = (mu-lambda*kbar)dt + sqrt(V)dZ + k dq
   dV = (alpha-beta*V)dt + sigmav *sqrt(v) d Zv

 cov(dZ,dZv)=corr(dZ,dZv)=rho dt

 prob(dq=1) = lambda dt,   ln(1+k) ~ N(ln(1+kbar) - 0.5 delta^2,delta^2)
 ___________________________________________________________________
 The risk neutral process is given as follows:

 dS/S = (b-lambdastar*kbarstar) dt + sqrt(V) dZstar + kstar dqstar
   dV = (alpha-betastar*V) dt + sigmav *sqrt(v) d Zv

 cov(dZ,dZv)=corr(dZ,dZv)=rho dt

 prob(dqstar=1) = lambdastar dt,   E(kstar)=kbarstar, var(kstar)=delta
 ----------------------------------------------------------------------
                      

 Note: values of b represent annualized rates for cost-to-carry
       - for non-div. paying stocks, indices or financial futures this is
         simply the rfr
       - for div paying stock or indices this is rfr-divrate
       - for currency options this is rfr(domestic)-rfr(foreign)

 This procedure must be accompanied by genoptfit.dec.
 Both procedures should be placed in your /gauss/src/ directory

*/
proc(1)=genoptfit(s0,v0,b,coc,tau,strike,theta);
/* note b refers to rfr, coc refers to rfr-divrate */
local maxrow,cvec,i,p1,p2;
local lowbd,upbd;
#include genoptfit.dec;
_intrec=0;
/* expand scalers to maximum vector length */
maxrow=maxc(rows(s0)|rows(v0)|rows(b)|rows(tau)|rows(strike));
if rows(s0)==1;
   s0=s0.*ones(maxrow,1);
endif;
if rows(v0)==1;
   v0=v0.*ones(maxrow,1);
endif;
if rows(b)==1;
   b=b.*ones(maxrow,1);
endif;
if rows(coc)==1;
   coc=coc.*ones(maxrow,1);
endif;
if rows(tau)==1;
   tau=tau.*ones(maxrow,1);
endif;
if rows(strike)==1;
   strike=strike.*ones(maxrow,1);
endif;
/* initialize vector of call option values */
cvec=zeros(maxrow,1);
_p1=cvec; _p2=cvec;
 /* set global parameters to pass into integration routines */
 _lambdastar=theta[1];
 _kbarstar=theta[2];
 _delta=theta[3];
 _alpha=theta[4];
 _betastar=theta[5];
 _sigv=theta[6];
 _rho=theta[7];
/* obtain call prices 1 at a time */
lowbd=0|10|100|1000; 
upbd=10|100|1000|10000; 
i=1;
do while i le maxrow;
 /* initialize other globals to pass into intergration routines */
 _v0=v0[i];
 _b=b[i];
 _tau=tau[i];
 _i=sqrt(-1);
 _s0=s0[i];
 _x=ln(strike[i]./s0[i]);
 /* obtain "probabilities" from moment generating functions */
 _mui=0.5;
 _intord=40;
 p1=0.5+(1./(pi)).*real(sumc(intquad1(&mgf,upbd'|lowbd')));
 _mui=-0.5;
 _intord=40;
 p2=0.5+(1./(pi)).*real(sumc(intquad1(&mgf,upbd'|lowbd')));
 _intord=40;
 cvec[i]=exp((coc[i]-_b).*_tau).*s0[i].*p1-exp(-_b.*_tau).*strike[i].*p2;
 _p1[i]=p1; _p2[i]=p2;
i=i+1;
endo;
retp(cvec); /*return vector of call prices */
endp;

proc(1)=mgf(phi);
local c,d,f;
phi=phi.*_i;
{c,d}=getcd(phi);
f=exp(c+d.*_v0+_lambdastar.*_tau.*((1+_kbarstar).^(_mui+0.5)).*
      (((1+_kbarstar).^phi).*exp((_delta.^2).*(_mui.*phi+0.5.*(phi.^2)))-1));
retp(((imag(f.*exp(-phi.*_x))))./(phi./_i));
endp;


proc(2)=getcd(phi);
local betai,rcoef2,coef1,coef0,coefc,c,d;
betai=_betastar-(0.5+_mui).*_rho.*_sigv;
rcoef2=0.5.*(_sigv.^2);
coef1=_rho.*_sigv.*phi-betai;
coef0=_mui.*phi+0.5.*(phi.^2);
coefc=(_b-_lambdastar.*_kbarstar).*phi;
disable;
if rcoef2 le (10.^-10);
 {c,d}=case1(_tau,coefc,_alpha,rcoef2,coef1,coef0);
else;
 {c,d}=case2(_tau,coefc,_alpha,rcoef2,coef1,coef0);
endif;
enable;
retp(c,d);
endp;

proc(2)=case1(t,coefc,alpha,rcoef2,coef1,coef0);
local c,d;
d=(coef1.==0).*(coef0.*t)+(coef1.ne 0).*(coef0./coef1).*(exp(coef1.*t)-1);
c=(coef1.==0).*(coefc.*t+0.5.*alpha.*coef0.*(t.^2))+(coef1.ne 0).*
  ((alpha.*coef0./(coef1.^2)).*(exp(coef1.*t)-1)-alpha.*(coef0./coef1).*t+
  coefc.*t);
retp(c,d);
endp;

proc(2)=case2(t,coefc,alpha,rcoef2,coef1,coef0);
local delta,gam,zdum,ratio3,egamt,cdenom,ddenom,d,c3dum,c3,c4,c;
delta=(coef1.^2)-4.*rcoef2.*coef0;
gam=sqrt(delta);
zdum=gam.<(10^-8);
ratio3=rcoef2.*coef0./(coef1.^2);
egamt=exp(gam.*t);
ddenom=coef1+(gam./(1-egamt)).*(1+egamt);
d=(zdum).*(2.*coef0.*t./(2-coef1.*t))
  +(1-zdum).*(-2.*coef0./ddenom);
c3dum=(ratio3.<(5*10^-6)).*(coef1.<0);
cdenom=rcoef2.*(1+ratio3+2.*(ratio3.^2));
c3=c3dum.*(coefc-alpha.*coef1./cdenom).*t+
   (1-c3dum).*(coefc-2.*alpha.*coef0./(coef1+gam)).*t;
c4=ln(1+0.5.*(1-coef1./gam).*(egamt-1));
c=(zdum).*((coefc-0.5.*alpha.*coef1./rcoef2).*t-(alpha./rcoef2).*ln(1-0.5.*
  coef1.*t))
 +(1-zdum).*(c3-(alpha./rcoef2).*c4);
retp(c,d);
endp;