Pricing Asian options with Monte Carlo Methods

[Robert Cyr]

With rather simple techniques, J can handle this problem with relative ease,
at least in one type of simulation. There is ample literature on this
subject, so I am only evoking the subject, with i hope a little J twist.
The calculations produce results similar to that obtained from the
followings source, so I assume they are correct.  I am submitting this
because it is fun(to me at least) and to find out if my technique can be
improved upon.

Referring again to Pierre l'Ecuyer's examples, who in turn refers us to a
book by P. Glasserman: Monte Carlo Methods in Financial Engineering.
Springer-Verlag,  2004, we find interesting assumption for the future
pricing of an asset.  This asset value simulation is the heart of the
problem.

First let's define succinctly the price of an asian option. In pricing an
asian option you are discounting at a risk-free rate the value of an asset
in excess of a certain amount(stk , the strike value).  With an asian
option, the value at a prescribed date is the average of the market value of
the security over an agreed period.

Then , the interesting part,  comes the simulation of the change in asset
value.  In this simulation, the projected asset value moves in accordance
with a geometric Brownian motion. Basically the log of the increase in value
is a linear function of time plus a normally distributed variation.

In this case the technique is quite transparently be applied to a 100,000
sample. This being only a demonstration the time periods considered appear
rather limited.

asian=: 3 : 0
'r sig stk s0 s n'=.0.05 0.5 100 100 12 100000
NB. r: risk free appreciation
NB. sig: volatility parameter
NB. stk: strike value
NB. s0: initial value
NB. s: time interval
NB. n: the sample size
dstat b=:(^-r)*stk-~stk>.(%s)*+/^}.+/\(^.s0),(s%~r--:*:sig)+(sig*%:%s)*normalrand
s,n
)

where
+/\(^.s0),(s%~r--:*:sig)+(sig*%:%s)*normalrand : the asset increases
exponentially, so its log is a linear function of time plus a normaly
distributed variation.
(%s)*+/^ is the average for calculation, the exponent changes from log to
value.
stk-~stk>. is the excess of this mean over the strike price(the agreed
future value)
(^-r)* is the discount factor


The timing of this routine is ok for one asset at 0.82 secs, but it adds up
and another technique is required where the number of cases is substantial.
----------------------------------------------------------------------
For information about J forums see http://www.jsoftware.com/forums.htm

[Mike Day-3]

It might help to analyse the routine with the timing monitor, JMP:
(it will also help to display this with a fixed width font)

   load'jpm'
   start_jpm_''
357142
   ts'asian 1'  NB. could just do bare call on "asian 1"
0.602233 1088
   showdetail_jpm_'asian'
 Time (seconds)
+--------+--------+---+--------------------------------------------------+
|all     |here    |rep|asian                                             |
+--------+--------+---+--------------------------------------------------+
|0.000030|0.000030|2  |monad                                             |
|0.000047|0.000047|2  |[0] 'r sig stk s0 s n'=.0.05 0.5 100 100 12 100000|
|0.817595|0.247245|2  |[7] z=.(sig%%:s)*normalrand s,n                   |
|0.029232|0.029232|2  |[8] z=.(s%~r--:*:sig)+z                           |
|0.087401|0.087401|2  |[9] z=.}.+/\(^.s0),z                              |
|0.145593|0.145593|2  |[10] z=.^z                                        |
|0.011769|0.011769|2  |[11] z=.(%s)*+/z                                  |
|0.007687|0.007687|2  |[12] b=:(^-r)*stk-~stk>.z                         |
|0.171057|0.000527|2  |[13] dstat b                                      |
|1.270410|0.529529|2  |total monad                                       |
+--------+--------+---+--------------------------------------------------+
   showdetail_jpm_'normalrand'
 Time (seconds)
+--------+--------+---+--------------------------------------+
|all     |here    |rep|normalrand                            |
+--------+--------+---+--------------------------------------+
|0.000016|0.000016|2  |monad                                 |
|0.790204|0.570867|2  |[0] (2 o.+:o.rand01 y)*%:-+:^.rand01 y|
|0.790220|0.570882|2  |total monad                           |
+--------+--------+---+--------------------------------------+
   showdetail_jpm_'rand01'
 Time (seconds)
+--------+--------+---+------------+
|all     |here    |rep|rand01      |
+--------+--------+---+------------+
|0.219338|0.219338|4  |monad ?@$ 0:|
+--------+--------+---+------------+
 
Evidently normalrand rand0 and dstat (!) are relatively heavy in
timing. There's been a lot of chat as you know about sampling
normal variates;  the cumulative sum and exponentiation of that
sum are also quite expensive.

Also
   showdetail_jpm_'dstat'
 Time (seconds)
+--------+--------+---+-----------------------------------------------------+
|all     |here    
|rep|dstat                                                |
+--------+--------+---+-----------------------------------------------------+
|0.000021|0.000021|2  
|monad                                                |
|0.000006|0.000006|2  |[0] t=.'/sample
size/minimum/maximum/median/mean'    |
|0.000010|0.000010|2  |[1] t=.t,'/std
devn/skewness/kurtosis'               |
|0.000032|0.000032|2  |[2] t=.,&':  ';._1
t                                 |
|0.000031|0.000031|2  |[3]
v=.$,min,max,median,mean,stddev,skewness,kurtosis|
|0.170442|0.000416|2  |[4] t,.":,.v
y                                       |
|0.170543|0.000517|2  |total
monad                                          |
+--------+--------+---+-----------------------------------------------------+
 
... suggests it might be worth seeing what makes this relatively slow!
Stick to min max mean & stddev?

Good luck!

Mike

Robert Cyr wrote:

> With rather simple techniques, J can handle this problem with relative ease,
> at least in one type of simulation. There is ample literature on this
> subject, so I am only evoking the subject, with i hope a little J twist.
> The calculations produce results similar to that obtained from the
> followings source, so I assume they are correct.  I am submitting this
> because it is fun(to me at least) and to find out if my technique can be
> improved upon.
>
> Referring again to Pierre l'Ecuyer's examples, who in turn refers us to a
> book by P. Glasserman: Monte Carlo Methods in Financial Engineering.
> Springer-Verlag,  2004, we find interesting assumption for the future
> pricing of an asset.  This asset value simulation is the heart of the
> problem.
>
> First let's define succinctly the price of an asian option. In pricing an
> asian option you are discounting at a risk-free rate the value of an asset
> in excess of a certain amount(stk , the strike value).  With an asian
> option, the value at a prescribed date is the average of the market value of
> the security over an agreed period.
>
> Then , the interesting part,  comes the simulation of the change in asset
> value.  In this simulation, the projected asset value moves in accordance
> with a geometric Brownian motion. Basically the log of the increase in value
> is a linear function of time plus a normally distributed variation.
>
> In this case the technique is quite transparently be applied to a 100,000
> sample. This being only a demonstration the time periods considered appear
> rather limited.
>
> asian=: 3 : 0
> 'r sig stk s0 s n'=.0.05 0.5 100 100 12 100000
> NB. r: risk free appreciation
> NB. sig: volatility parameter
> NB. stk: strike value
> NB. s0: initial value
> NB. s: time interval
> NB. n: the sample size
> dstat b=:(^-r)*stk-~stk>.(%s)*+/^}.+/\(^.s0),(s%~r--:*:sig)+(sig*%:%s)*normalrand
> s,n
> )
>
> where
> +/\(^.s0),(s%~r--:*:sig)+(sig*%:%s)*normalrand : the asset increases
> exponentially, so its log is a linear function of time plus a normaly
> distributed variation.
> (%s)*+/^ is the average for calculation, the exponent changes from log to
> value.
> stk-~stk>. is the excess of this mean over the strike price(the agreed
> future value)
> (^-r)* is the discount factor
>
>
> The timing of this routine is ok for one asset at 0.82 secs, but it adds up
> and another technique is required where the number of cases is substantial.
> ----------------------------------------------------------------------
> For information about J forums see http://www.jsoftware.com/forums.htm
>
> No virus found in this incoming message.
> Checked by AVG - http://www.avg.com 
> Version: 8.0.138 / Virus Database: 270.6.3/1613 - Release Date: 15/08/2008 05:58
>
>
>
> No virus found in this incoming message.
> Checked by AVG - http://www.avg.com 
> Version: 8.0.138 / Virus Database: 270.6.3/1613 - Release Date: 15/08/2008 05:58
>
>
>
>  

----------------------------------------------------------------------
For information about J forums see http://www.jsoftware.com/forums.htm