@===========================================================================@
@ Option Pricing with Stochastic Volatility: Hull-White Model with rho=-0.5 @
@===========================================================================@

new;
S=100;           @ spot price @
X=100;           @ exercise or strike price @
tau=0.5;         @ time to maturity @
r=0.1;           @ risk-free rate @
n=90;            @ number of subperiods @
nmc=5000;        @ number of Monte Carlo simulations @
dt=tau/n;        @ time step @
sigv=1;          @ volatility of volatility @
a=10;            @ speed of mean reversion @
sigstar=(.15)^2; @ long-run mean level of volatility @
V0=(.1)^2;       @ initial value for volatility @
rho=-0.5;        @ correlation between stock price and volatility @
bsc=bs_call(S,X,tau,r,sqrt(V0));

@ Monte Carlo Simulation with Antithetic and Control Variates @
callsv=zeros(nmc,1); p1=zeros(nmc,1); p2=zeros(nmc,1);
p3=zeros(nmc,1); p4=zeros(nmc,1); q1=zeros(nmc,1); q2=zeros(nmc,1);
i=1; do while i<=nmc;
 V1=V0*ones(n+1,1);  V2=V0*ones(n+1,1);
 S1=S*ones(n+1,1);  S2=S*ones(n+1,1); S3=S*ones(n+1,1);
 S4=S*ones(n+1,1); S5=S*ones(n+1,1);  S6=S*ones(n+1,1);
 sigma=(1~rho)|(rho~1);
 err=rndn(n+1,2)*chol(sigma); err2=-err;
 j=2; do while j<=n+1;
  V1[j]=V1[j-1]*exp((a*(sigstar-V1[j-1])-(sigv^2)/2)*dt+err[j,1]*sigv*sqrt(dt));
  V2[j]=V2[j-1]*exp((a*(sigstar-V2[j-1])-(sigv^2)/2)*dt+err2[j,1]*sigv*sqrt(dt));
  S1[j]=S1[j-1]*exp((r-V1[j]/2)*dt+err[j,2]*sqrt(V1[j]*dt));
  S2[j]=S2[j-1]*exp((r-V1[j]/2)*dt+err2[j,2]*sqrt(V1[j]*dt));
  S3[j]=S3[j-1]*exp((r-V2[j]/2)*dt+err[j,2]*sqrt(V2[j]*dt));
  S4[j]=S4[j-1]*exp((r-V2[j]/2)*dt+err2[j,2]*sqrt(V2[j]*dt));
  S5[j]=S5[j-1]*exp((r-V0/2)*dt+err[j,2]*sqrt(V0*dt));
  S6[j]=S6[j-1]*exp((r-V0/2)*dt+err2[j,2]*sqrt(V0*dt));
 j=j+1; endo;
p1[i]=exp(-r*tau)*maxc((S1[rows(S1)]-X)|0);
p2[i]=exp(-r*tau)*maxc((S2[rows(S2)]-X)|0);
p3[i]=exp(-r*tau)*maxc((S3[rows(S3)]-X)|0);
p4[i]=exp(-r*tau)*maxc((S4[rows(S4)]-X)|0);
q1[i]=exp(-r*tau)*maxc((S5[rows(S5)]-X)|0);
q2[i]=exp(-r*tau)*maxc((S6[rows(S6)]-X)|0);
callsv[i]=((p1[i]+p3[i]-2*q1[i])/2+(p2[i]+p4[i]-2*q2[i])/2)/2;
i=i+1; endo;

? "Price of a call option with stochastic vol. ";;bsc+meanc(callsv);
? "BS price of a call option with constant vol.";;bsc;

@ Black-Scholes Formula for a Call Option @
proc bs_call(S,X,T,r,sigma);
  local d1,d2,c;
  d1=(ln(S/X)+(r+(sigma^2)/2)*T)/(sigma.*sqrt(T));
  d2=d1-sigma.*sqrt(T);
  c=S.*cdfn(d1)-X.*exp(-r.*T).*cdfn(d2);
  retp(c);
endp;