(* :Title: EulerSimulate *)

(* :Author: Mark Fisher *)

(* :Context: EulerSimulate` *)

(* :Mathematica Version: 4.0 *)

(* :History: 
		Version 1.0 February 1999.
		Version 1.1 September 1999.
			Use RandomArray instead of StandardNormal.
		Version 1.2 June 2000.
			Improved compiling.
			Modified argument list.
			Added option InitialTime.
			Added scalar verion.
*)

(* :Keywords: Ito's lemma, stochastic calculus, random walk, brownian
motion, option pricing, Black-Scholes *)

(* :Summary: simulate SDEs via Euler approximation *)

(* :Discussion:
This package impliments first-order Euler simulation for any number of Ito
processes with any number of arbitrarily-correlated Brownian motions. The
sole function is EulerSimulate.
*)

(* :Requirements: This package imports the package 
	Statistics`NormalDistribution.m *)

(* :Examples:

Example 1
Geometric brownian motion:

sim = EulerSimulate[.2 s, .1 s, {s, 100}, {1, 10^3}];
ListPlot[sim, PlotJoined -> True]

Example 2
To let time play an active role, model time explicitly as an Ito process:

sim = EulerSimulate[{1, x - t}, {{0}, {1}}, {t, 0}, {x, 0}, {1, 10^3}, 
	IncludeTime -> False];
ListPlot[sim, PlotJoined -> True]

Example 3
Three "affine" state variables:

mu = {5(m - r), .1( .05 - m), 1( .01 - v)};
sig = Re @ {{Sqrt[v], 0, 0}, {0, .1 Sqrt[m], 0}, {0, 0, .1 Sqrt[v]}};
sim = EulerSimulate[mu, sig, {r, .05}, {m, .05}, {v, .01}, {1, 10^3}];
Show @ GraphicsArray @
	Table[ListPlot[sim[[All, {1, i}]], PlotJoined -> True,
		DisplayFunction -> Identity, Frame -> True, Axes -> False], 
		{i, 2, 4}]

*)

BeginPackage["EulerSimulate`", {"Statistics`NormalDistribution`",
	"Statistics`Common`DistributionsCommon`"}]

EulerSimulate::usage = "EulerSimulate[drift, diffusion, {x, x0}, {duration,
nsteps}] returns a list of simulated values for the Ito process determined
by the arguments. The arguments specify the drift and diffusion, the
symbolic name and initial value {x, x0}, and the length of time and the
number of time steps,{duration, nsteps}. The default settings for the
options are Compiled -> True, IncludeTime -> True, InitialTime -> 0, and
IntermediateValues -> True. A system of Ito processes can be simulated by
specifying the drift as a vector and the the diffusion as a matrix, and by
providing a {x, x0} pair for each process in a sequence. By default,
EulerSimulate generates a matrix of orthogonal standard normal deviates;
alternatively, the vector form of EulerSimulate allows the user to supply a
matrix of shocks as a final optional argument."

IncludeTime::usage = "IncludeTime is an option for EulerSimulate. The
default setting is IncludeTime -> True."

InitialTime::usage = "InitialTime is an option for Eulersimulate. The
default setting is InitialTime -> 0."

IntermediateValues::usage = "IntermediateValues is an option for
EulerSimulate. The default setting is IntermediateValues -> True."

Begin["`Private`"]

Options[EulerSimulate] = {Compiled -> True, IncludeTime -> True,
	InitialTime -> 0, IntermediateValues -> True}

(* main engine: vector version, takes shocks as argument *)

EulerSimulate[drift_?VectorQ, diffusion_?MatrixQ,
	x:{_Symbol, _?NumericQ}..,
	{duration_?NumericQ, nsteps_Integer},
	devs_?(MatrixQ[#, NumberQ]&), opts___?OptionQ] /;
		(Length[drift] == Length[diffusion] == Length[{x}] &&
		Dimensions[devs][[2]] == Dimensions[diffusion][[2]]) :=
	Module[{compile, time, t0, list, dt, sqrtdt, args0, args, 
		errargs, fun, sim},
	{compile, time, t0, list} = 
		{Compiled, IncludeTime, InitialTime, IntermediateValues} /. 
		{opts} /. Options[EulerSimulate];
	sqrtdt = Sqrt[ dt = N[duration/nsteps] ];
	{args, args0} = Transpose[{x}];
	errargs = Unique[]& /@ Range[Last @ Dimensions[diffusion]];
	fun = If[TrueQ[compile], Compile, Function] @@ 
		{Join[args, errargs], args + drift dt + (diffusion.errargs) sqrtdt};
	sim = If[TrueQ[list], FoldList, Fold]
		[fun @@ Join[#1, #2]&, N[args0], devs];
	If[TrueQ[time] && TrueQ[list],
		Flatten /@ Transpose[
			{Range[N[t0], t0 + duration, duration/nsteps], sim}],
		sim
		]
	]

(* vector version, constructs shocks and calls main engine *)

EulerSimulate[drift_?VectorQ, diffusion_?MatrixQ,
	x:{_Symbol, _?NumericQ}..,
	{duration_?NumericQ, nsteps_Integer}, opts___?OptionQ] /;
	(Length[drift] == Length[diffusion] == Length[{x}]) :=
		EulerSimulate[drift, diffusion, x, {duration, nsteps},
			(* make the random number matrix *)
			RandomArray[NormalDistribution[0,1], 
				{nsteps, Last @ Dimensions[diffusion]}], 
			opts]

(* scalar version, calls vector version *)

EulerSimulate[drift_, diffusion_, {x_Symbol, x0_?NumericQ}, 
		{duration_?NumericQ, nsteps_Integer}, opts___?OptionQ] :=
	If[TrueQ[IncludeTime /. {opts} /. Options[EulerSimulate]],
		(* then do nothing *)
		Identity, 
		(* else flatten the result *)
		Flatten] @ 
		EulerSimulate[{drift}, 
			If[VectorQ[diffusion], {diffusion}, {{diffusion}}], 
			{x, x0}, {duration, nsteps}, opts]

End[]
EndPackage[]