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

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 @

@ Monte Carlo Simulation with antithetic variates @
callp=zeros(nmc,1);
i=1; do while i<=nmc;
 v1=V0*ones(n+1,1); v2=V0*ones(n+1,1);
 j=2; do while j<=n+1;
  nu1=rndn(1,1); nu2=-nu1;
  v1[j]=v1[j-1]*exp((a*(sigstar-v1[j-1])-(sigv^2)/2)*dt+nu1*sigv*sqrt(dt));
  v2[j]=v2[j-1]*exp((a*(sigstar-v2[j-1])-(sigv^2)/2)*dt+nu2*sigv*sqrt(dt));
 j=j+1; endo;
callp1=bs_call(S,X,tau,r,sqrt(meanc(v1[2:n])));
callp2=bs_call(S,X,tau,r,sqrt(meanc(v2[2:n])));
callp[i]=(callp1+callp2)/2;
i=i+1; endo;

? "Price of a call option with stoch.volatility";;meanc(callp);
? "Black-Scholes price with constant volatility";;bs_call(S,X,tau,r,sqrt(V0));

@ 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;