@ File STABOPT @
@ J. Huston McCulloch 4/15/03 @ 
@ Computes Stable Option Values From Risk-Neutral Measure 
   by method of "The Risk-Neutral Measure under Log-Stable 
   Uncertainty," econ.ohio-state.edu/jhm/papers/rnm.pdf @ 
@ Contains procs 
    STABOPT -- computes stable option values by Fourier inversion
    TAILOPT -- computes far OTM stable option values
    STABPHI -- computes Phi function for TAILOPT
    INCGAMMA -- computes incomplete gamma fn for any a, except nonpositive integers
    STABRNMCF -- computes CF of Stable RNM for ROMBFFTI  
    ROMBFFTI -- gives centered Romberg FFTI of centered function   
    FFTROT -- Gives centered FFT of centered function (uses GAUSS "FFTI")
    FFTIROT -- Gives centered FFTI of centered function (uses GAUSS "FFT")    
    (ROMBFFT not implemented, but straightforward) 
    SGN -- computes sign of real part of complex t @
    
@ Illustrates STABOPT @
new;
library pgraph;
graphset;
_pdate = "";
_pltype = 6;
_plwidth = {12,6};
_pgrid = {1,0};
_pcolor = {1, 12, 2, 5}; @ = dk blue, light red, dk green, dk magenta @ 
_pnum = 2; @ turns y axis numbers @
_pmcolor = 0; @ sets axes black @
_pnumht = .20;
_paxht = .20; 
_ptitlht = .20;
fonts("microb");
format /rdn 7,3;

s0 = 100;       @ spot price @ 
r = .05;        @ annual interest rate @
div = .03;      @ annual dividend rate @ 
t = 1;          @ years to maturity @
alpha = 1.5;
beta =0;
c0 = .1;        @ scale that accumulates in 1 year @
n = 2^10;       @ must be a power of 2 @
nx = 200;
z = seqa(ln(.5), ln(2)/(nx/2),nx);
x = exp(z)*s0; @ vector of nx exercise prices @ 

{cv, pv, e} = stabopt(s0,x,r,div,t,alpha,beta,c0,n);
votm = minc((cv~pv)'); @ computes OTM option values from 
  call values CV, put values PV @
title ("OTM option values, alpha = 1.5, beta = 0");
xlabel("Exercise Price");
ylabel("Option Value");
xy(x,votm~(votm+e)); @ second line differs by estimated measurement error @ 
"waiting"; 
wait;

@ Compute and plot Volatility Smile @ 
smile = zeros(nx,1);
for i (1,nx,1);
 smile[i] = EuropeanBSCall_ImpVol(cv[i],s0,x[i],r,div,t);
 endfor;
ytics(0,.5,.1,5);
title ("Volatility Smile, alpha= 1.5, beta = 0");
xlabel ("Log Exercise Price");
ylabel ("BS Implied Volatility (s.d.)");
xy(z,smile);
end;

proc(3) = stabopt(s0, x, r, div, t, alpha, beta, c0, n);
@ computes stable option prices using Romberg FFT inversion@
@ S0 scalar is current asset price
  x is column vector of increasing, nonnegative exercise prices
  r scalar is risk-free annual (or periodic) 
   interest rate (fraction, not %)
  div scalar is annual (or periodic) dividend payout rate 
   (fraction, not %)
  t scalar is remaining time to maturity in years (or periods) 
  alpha scalar must be in [1.3, 2]
  beta scalar must be in [-1,1]
  c0 scalar is stable scale that accumulates in 1 year (or 1 period) 
  n scalar is # of points for FFT.  Must be even power of 2. 
   n = 2^10 or less is fast, but check errors.  
   n = 2^14 may take 10 seconds, but gives good precision @
@ Returns vc, vp, e 
  vc = vector value of calls 
  vp = vector value of puts
  abs(e) is vector estimate of the FFT integration error.  
   e > 0 indicates FFT values are used.
   e < 0 indicates tail value has been substituted for 
    FFT value.  This occurs when difference is less than 
    abs(e). 
   e = 0 indicates FFT is out of range and tail value 
    is used. @  
if alpha < 1.3;
  print "Alpha must be >= 1.3 in STABOPT for adequate precision";
  print "Actual alpha value is "; alpha;
  stop;
  endif;
if not x==sortc(x,1);
  print "X must be in increasing order for STABOPT";
  stop;
  endif;    
local z, vz, ez, xz, mask, ind, th, v, e, vt, tail;
local zx, dz, cz, inrange, outrange, nzx, vinrange;
local vc, vp, einrange, a1, a2, a3, f, c; 
f = s0*exp((r-div)*t);
x = x/f;
c = c0*t^(1/alpha);
declare astabopt_, bstabopt_, cstabopt_;
astabopt_ = alpha;
bstabopt_ = beta;
cstabopt_ = c;
{z,vz,ez} = ROMBFFTI(&stabphi,n,c);
vz = abs(vz.*(real(vz).>0)); 
xz = exp(z);
cz = vz.*(xz.>=1) + (vz+1-xz).*(xz.<1);
 @ gives call values at all n z values @
inrange = (x.>xz[3]) .and (x.<xz[n-2]);
outrange = indexcat(inrange,0);
inrange = indexcat(inrange,1);
zx = ln(x[inrange]);
nzx = rows(zx);
ind = sortc(sortind(zx|z)~seqa(1,1,nzx+n),1);
ind = ind[1:nzx,2] - seqa(1,1,nzx);
@ ind now contains the indices of the z values immediately 
  below the zx values @
dz = z[2]-z[1];
@ for internal x's, @ 
th = (zx-z[ind])/dz;
a2 = (cz[ind-1] + cz[ind+1])/2 - cz[ind];
a3 = (-cz[ind-1] + 3*cz[ind] - 3*cz[ind+1] + cz[ind+2])/6;
a1 = (cz[ind+1]-cz[ind-1])/2 - a3;
vc = cz[ind] + th.*(a1 + th.*(a2 + th.*a3));
vinrange = vc.*(zx.>=0) + (vc + x[inrange] -1).*(zx.<0); 
einrange = ez[ind].*(1-th) + ez[ind+1].*th;
vt = tailopt(alpha, beta, c, x);
tail = abs(vt[inrange]-vinrange).<einrange .and vt[inrange] ./= -1;
v = zeros(rows(x),1);
v[outrange] = vt[outrange];
v[inrange] = vinrange.*(.not tail) + vt[inrange].*tail;
e = zeros(rows(x),1);
e[inrange] = einrange.*(.not tail - tail);
vp = v.*(x.<1) + (v+x-1).*(x.>=1);
vc = (v+1-x).*(x.<1) + v.*(x.>=1);
vp = f * vp * exp(-r*t);
vc = f * vc * exp(-r*t);
e = f * e * exp(-r*t);
retp(vc, vp, e);
endp;

/*
@ tests stabphi @
astabopt_ = .9;   
  @ alpha = .9 gives wrong sign on real part with beta = 1 @
  @ present application requires alpha >= 1.3, so moot for now @            
bstabopt_ = 1;
cstabopt_ = 1;
print astabopt_ bstabopt_ cstabopt_;
t = seqa(-1,.001,2000);
phi = stabphi(t);
xtics(-1,1,.5,5);
xy(t,abs(phi)~real(phi)~imag(phi));
end;
*/

proc stabphi(t);
@ computes stable OTM option value phi function @
@ receives alpha, c, beta as globals astabopt_, etc. @
@ assumes F = 1 @
@ not implemented for alpha = 1 @
local aye, phi, ind, theta, alpha, c, beta;
aye = sqrt(-1);
alpha = astabopt_;
beta = bstabopt_;
c = cstabopt_;
if abs(alpha-1)<1.e-8;
 print "Stabphi not implemented for alpha = 1";
 stop;
 endif;
theta = pi*alpha/2;
ind = abs(t).<1.e-8; 
phi = (stabrnmcf(t-aye)-1)./(aye*(t+ind).*(aye*t+1));
phi = phi.*(1-ind)-ind.*(c^alpha)*(beta + .5*(1-beta)*alpha)/cos(theta);
retp(phi);
endp;

proc stabrnmcf(t); 
@ computes stable RNM CF for STABOPT @ 
@ receives alpha, c, beta as globals astabopt_, etc.  @ 
@ assumes F = 1 @
local aye, rnm, theta, alpha, c, beta;
aye = sqrt(-1);
alpha = astabopt_;
beta = bstabopt_;
c = cstabopt_;
if abs(alpha-1)<1.e-8;
 print "Stabrnmcf not implemented for alpha = 1";
 stop;
 endif;
theta = pi*alpha/2;
rnm = -(aye*beta/cos(theta))*t ;
rnm = rnm - (.5*(1-beta)*(t.*sgn(t)).^alpha).*(1+aye*sgn(t)*tan(theta)); 
rnm = rnm + .5*(1+beta).*(1-(1-aye*t).^alpha)/cos(theta);
rnm = exp((c^alpha)*rnm);
retp(rnm);
endp;

proc tailopt(alpha, beta, c, x);
local n, z, v, psi, calls, puts, k, th;
n = rows(x);
if abs(alpha-1) < 1.e-5;
 th = (alpha - 1)/(2.e-5);  
 v = tailopt(1-1.e-5,beta,c,x)+th.*tailopt(1+1.e-5,beta,c,x);
 retp(v);
elseif alpha == 2;
 v = zeros(n,1);
 retp(v); 
 endif;
@ else @
v = -ones(n,1);
z = ln(x);
psi = v;
calls = indexcat(x.>1,1);
puts = indexcat(x.<1 .and x.>0 ,1);
k = gamma(alpha)*sin(pi*alpha/2)/pi;  
if calls /= miss(1,1);
psi[calls] = k*x[calls].*incgamma(1-alpha,z[calls]);
v[calls] = (c^alpha)*(1+beta)*psi[calls];
endif;
if puts /= miss(1,1);
psi[puts] = (k./x[puts]).*incgamma(1-alpha,-z[puts]);
v[puts] = (c^alpha)*(1-beta)*x[puts].*psi[puts];
endif;
v = v.*(x./=0);
retp(v);
endp;

proc incgamma(a,x);
@ Incomplete gamma function for a not a 
   nonpositive integer, by repeated 
   recursions if necessary for a < 0 @
if a > 0;
retp(gamma(a)*(1-cdfgam(a,x)));
else; 
retp((incgamma(a+1,x)-(x.^a).*exp(-x))/a);
endif;
endp;
 
@ Illustrates ROMBFFTI with recovery of Laplace 
  density from its cf @ 
new;
fn cf(t) = 1/(1+t.^2);
{x,y,e} = ROMBFFTI(&cf,2^10,1);
den = .5*exp(-abs(x));
print maxc(abs(y-den));
end;

proc(3) = ROMBFFTI(&func,n,c);
 @ Performs Romberg Inverse Fast Fourier Transform 
   of function func, using n, 4*n, and 16*n points, 
   and extrapolating to infinity.
   n must be an even power of 2.   
   c is a scale factor for the returned function 
     (eg the stable scale parameter).
   Computes approximate value of  
     y(x) = (1/pi) integral(-inf to inf) exp(-i x t) func(t) dt,
        where i = sqrt(-1).     
   Returns x, y, e.  
   x contains the n values at which inverse transform 
      y(x) is evaluated.  
   x is in steps of dx = c*sqrt(2*pi/n),  
      and takes values dx*(-n/2)... dx*(n/2 - 1).  
   y is not constrained to be real or nonnegative. 
     Set y = abs(y) if you are expecting a real return.  
   e is an estimate of error based on starting with n/4 points.     
     For the middle half of the returned x values, e is 
     equal to its value at an adjacent point.  For the outer half, 
     e equals its value at the outermost value in the inner half. 
   The integration of func(t) is performed in n steps of 
     dt = 2*pi/(n*dx) rads/xunit   @

local dx, dt, fcrit, xmin, x, t, f, y, y1, x1;
local y2, y3, d2, d3, ind, rho;
local y0, e, d1, n0, yy; 
local func:fn; 

if 2^trunc(ln(n+.5)/ln(2)) /= n;
 print "n must be an even power of 2 in ROMBFFTI";
 print "Actual value of n is " n;
 stop;
 endif;
if n < 4; 
 print "n must be > 4";
 print "Actual value of n is " n;
 stop;
 endif;

n=n/4; 
for i (0, 3, 1);
dx = c*sqrt(2*pi/n);
dt = 2*pi/(n*dx); @rads/xunit@
fcrit = pi/dx; @rads/xunit@
xmin = -dx * n / 2;
x = seqa(xmin,dx,n);
t = seqa(-fcrit,dt,n);
f = func(t);
y = fftroti(f,dx);
y = abs(y);
if i == 0;
 y0 = y;
elseif i == 1; 
 y1 = y;
 x1 = x;
elseif i == 2;
 y2 = y;
else; 
 y3 = y;
endif;
n = n*4;
endfor;

n = rows(y1);
y2 = reshape(y2,2*n,2);
y2 = y2[n/2+1:3*n/2,1];
y3 = reshape(y3,4*n,4);
y3 = y3[3*n/2+1:5*n/2.,1];
d3 = y3-y2;
d2 = y2-y1;
ind = (abs(d3).>1.e-8).*(abs(d2).>1.e-8);  
    @ indicates neither d3 nor d2 is small@
rho = ind.*d3./(d2.*ind + (1-ind)); 
    @ = d3/d2 if d2 and d3 large, 0 if either small already @
y = y3 + d3.*rho./(1-rho);

@ estimate error e by doing same procedure for n/4 @
n0 = n/4;
y1 = reshape(y1,2*n0,2);
y1 = y1[n0/2+1:3*n0/2,1];
y2 = reshape(y2,2*n0,2);
y2 = y2[n0/2+1:3*n0/2,1];
d2 = y2-y1;
d1 = y1-y0;
ind = (abs(d2).>1.e-8).*(abs(d1).>1.e-8);  
rho = ind.*d2./(d1.*ind + (1-ind)); 
y0 = y2 + d2.*rho./(1-rho);
yy = reshape(y,2*n0,2);
yy = yy[n0/2+1:3*n0/2,1];
e = abs(y0-yy);
e = reshape(e~e,2*n0,1);
e = ones(n0,1)*e[1]|e|ones(n0,1)*e[2*n0];
retp(x1,y,e);
endp;

proc fftrot(y,dx);
@ Computes FFT of y(x) for x = -(n/2)dx, ... (n/2 - 1)dx @
@ Returns z(t) for t = -fcrit, ... fcrit(n/2 - 1)/(n/2), fcrit = pi/dx @
@ Gauss "ffti" actually computes FFT & vice-versa @
local n, yy, zz, z;
n = rows(y);
yy = y[n/2+1:n]|y[1:n/2];
zz = ffti(yy)*dx; 
z = zz[n/2+1:n]|zz[1:n/2];
retp(z);
endp;

proc fftroti(z,dx);
@ Computes FFTI of z(t) for t = -fcrit, ... fcrit(n/2 - 1)/(n/2), fcrit = pi/dx @
@ Returns y(x) for x = -(n/2)dx, ... (n/2 - 1)dx @
@ Gauss "fft" actually computes FFTI & vice-versa @
local n, zz, yy, y;
n = rows(z);
zz = z[n/2+1:n]|z[1:n/2];
yy = fft(zz/dx); 
y = yy[n/2+1:n]|yy[1:n/2];
retp(y);
endp;

proc sgn(t);
@ computes sign of real part of t 
  returns 0 if real part is 0 @ 
local s;
s = (real(t).>0)-(real(t).<0);
retp(s);
endp;