(* ::Package:: *)

(************************************************************************)
(* This file was generated automatically by the Mathematica front end.  *)
(* It contains Initialization cells from a Notebook file, which         *)
(* typically will have the same name as this file except ending in      *)
(* ".nb" instead of ".m".                                               *)
(*                                                                      *)
(* This file is intended to be loaded into the Mathematica kernel using *)
(* the package loading commands Get or Needs.  Doing so is equivalent   *)
(* to using the Evaluate Initialization Cells menu command in the front *)
(* end.                                                                 *)
(*                                                                      *)
(* DO NOT EDIT THIS FILE.  This entire file is regenerated              *)
(* automatically each time the parent Notebook file is saved in the     *)
(* Mathematica front end.  Any changes you make to this file will be    *)
(* overwritten.                                                         *)
(************************************************************************)



Protect[t,Z];  (* Save these symbols... *)

WeinerProcessMake::usage="WeinerProcessMake[s_Symbol] builds a Weiner
	process diffusion.  The procedure has no explicit result.";

DiffusionWeinerProcessMake::usage="DiffusionWeinerProcessMake[s_Symbol]
	builds a Weiner	process diffusion. The procedure has no explicit
	result and is an alias for WeinerProcessMake.";
	
	
DiffusionMake::usage="DiffusionMake[x_Symbol,wp_diffusion,f,g,x0] builds
	a diffusion	based on the Weiner process defined by wp.  The resulting
	diffusion has symbolic name x and the drift and dispersion functions
	f and g, respectively.  x0 is the intended intial value.
	The procedure has no explicit result.";

diffusion::usage="Internal data structure for holding the contents
	of a diffusion object.";
	

DiffusionLabel::usage="DiffusionLabel[d_diffusion] returns the
	text name of the diffusion d.";
	
DiffusionSymbol::usage="DiffusionSymbol[d_diffusion] returns the
	symbolic name of the diffusion d wrapped in HoldForm to avoid
	evaluation recursion.";
	
DiffusionDrift::usage="DiffusionDrift[d_diffusion] returns the drift
	function of	the	diffusion d. DiffusionDrift[d_diffusion,Function]
	returns a function object.";

DiffusionDispersion::usage="DiffusionDispersion[d] returns the dispersion
	function of	the	diffusion d. DiffusionDispersion[d_diffusion,Function]
	returns a function object.";
	
DiffusionInitialValue::usage="DiffusionInitialValue[d_diffusion] returns
	the initial value of the diffusion d.";

DiffusionPrint::usage = "DiffusionPrint[d] prints a formatted version
	of the diffusion object d using subscripts and symbolic names.";
	
DiffusionExpand::usage = "DiffusionExpand[expr] attempts to combine
	the diffusion elements in the input expression into a single
	diffusion with simple algebraic operations.";
	
DiffusionSimulate::usage="DiffusionSimulate[d_diffusion, n, dt] simulates a
	diffusion.  It returns a list of n pairs {,,,{i*dt,x[i*dt]},,,} of the
	diffusion d. Each realized series is an independent realization.
	Responds with an error if the diffusion contains symbolic paramters.";

DiffusionA::usage = "DiffusionA[expr] is an alias for A[d][expr].";

A::usage = "A[expr] applies the infinitesimal operator
	associated with the diffusion d to the expression expr.";

DiffusionIto::usage = "DiffusionIto[expr] is an alias for
	Ito[expr].";
	
Ito::usage = "Ito[expr] finds the diffusion that results
	from the scalar expression which is assumed to contain
	a diffusion.  Options: 
			'itoInvert -> True' --> default
			'itoSymbol -> Y   ' --> symbol for new process.";

ItoInvert::usage="If true (default), attempts to symbolically find
	the inverse diffusion.";
	
ItoSymbol::usage="New symbol denoting the diffusion created by 
	applyinhg the Ito formula.";



Clear[findDiffusion, replaceDiffusion, replaceHeld];

	findDiffusion::tooFew =
		 "This expression does not contain a diffusion.";
	findDiffusion::tooMany =
		 "This expression has more than one diffusion.";

	findDiffusion[expr_] := Module	\
		[	{found, len},
			found = Cases[expr, _diffusion, Infinity];
			len = Length[found];
			If[ len == 0, Message[findDiffusion::tooFew]; Return];
			If[ len > 1 , Message[findDiffusion::tooMany]; Return];
			First[found]
		]
				
	replaceDiffusion[expr_] :=
			expr /. d_diffusion :> DiffusionSymbol[d]; 
			
	replaceDiffusion[expr_,d_diffusion] := 
			expr /. d -> DiffusionSymbol[d]
					 
	replaceDiffusion[expr_,d_diffusion,sub_] := Module	\
		[	{hf = DiffusionSymbol[d]},
			expr /. {d -> sub, hf -> sub}
		]

	replaceHeld[expr_,str_,sym_Symbol] := 
			expr /. (HoldForm @@ ToHeldExpression[str]) -> sym


Clear[DiffusionMake, WeinerProcessMake, DiffusionChangeSymbol];

WeinerProcessMake[var_Symbol] := Module \
	[ 	{},
		Set[var, diffusion[HoldForm[var],HoldForm[var],0,1,0]];
	]

DiffusionMake[var_Symbol,diffusion[wp_,wp_,___], f_,g_,initial_:0] := Module \
	[	{drift, disp, iv},
		drift = f /. var -> HoldForm[var];
		disp  = g /. var -> HoldForm[var];
		iv    = initial /. var -> HoldForm[var];
		Set[var, diffusion[HoldForm[var],wp,drift,disp,iv]];
	]
	
DiffusionChangeSymbol[diffusion[x_,wp_,f_,g_,i_], new_Symbol] :=
	Module[ {},
		DiffusionMake[new,ReleaseHold[wp],
					f/.x->new,g/.x->new,i/.x->new];
		new];


Clear[DiffusionQ, HasHeldDiffusionQ]

DiffusionQ[x_diffusion] := True
DiffusionQ[_] := False

HasHeldDiffusionQ[_Symbol] := False;
HasHeldDiffusionQ[_?NumberQ] := False;
HasHeldDiffusionQ[HoldForm[_?DiffusionQ]] := True
HasHeldDiffusionQ[expr_] := 
		And @@ HasHeldDiffusionQ /@
			Cases[expr,HoldForm[_],Infinity]
HasHeldDiffusionQ[_] := False


Clear[diffusion, diffusionNewName, DiffusionExpand]


diffusionNewName := HoldForm[Evaluate[Unique[]]]

DiffToSymbolsQ = True;

Attributes[DiffusionExpand] = {HoldAll};

DiffusionExpand[expr_] := Block	\
	[	{res},
		DiffToSymbolsQ = False;
		res = Release[expr] //. DiffusionExpandRules;
		DiffToSymbolsQ = True;
		res
	]

DiffusionExpandRules = 
	{ k_?NumberQ * diffusion[_,w_,mu_,sig_,x0_] -> 
			diffusion[diffusionNewName,w,k mu, k sig, k x0],
	  k_Symbol   * diffusion[_,w_,mu_,sig_,x0_] -> 
			diffusion[diffusionNewName,w,k mu, k sig, k x0],
	  diffusion[_,w_,m1_,s1_,x01_] + diffusion[_,w_,m2_,s2_,x02_] ->
	  		diffusion[diffusionNewName,w,m1+m2,s1+s2,x01+x02] };


diffusion /:  k_ * d_diffusion :=
					k * DiffusionSymbol[d] /; DiffToSymbolsQ;
					
diffusion /:  k_ + d_diffusion :=
					k + DiffusionSymbol[d] /; DiffToSymbolsQ;
					
diffusion /:  f_[d_diffusion,t_] := 
				DiffusionSymbol[f] @@ {DiffusionSymbol[d],t} /;
					(Head[f]===diffusion) && DiffToSymbolsQ;


Clear[DiffusionSymbol,
	  DiffusionWeinerProcess, DiffusionWeinerProcessSymbol,
	  DiffusionDrift,
	  DiffusionDispersion, DiffusionInitialValue];
	  
Begin["PrivateDiffAccess`"];


DiffusionWeinerProcess [diffusion[_,wp_,___]]  := ReleaseHold[wp]

DiffusionSymbol[diffusion[l_,___]] := l;
DiffusionWeinerProcessSymbol[diffusion[_,wp_,___]] :=
		DiffusionSymbol[ReleaseHold[wp]];


DiffusionDrift[diffusion[_,_,f_,___]]  := f;
DiffusionDrift[d_diffusion, Function] := 
	Function[{x,t}, Release[replaceDiffusion[DiffusionDrift[d],d,x]]  ]

DiffusionDispersion[diffusion[_,_,_,s_,___]]  := s;
DiffusionDispersion[d_diffusion, Function] := 
	Function[{x,t}, Release[replaceDiffusion[DiffusionDispersion[d],d,x]]  ]
	
DiffusionInitialValue[diffusion[_,_,_,_,iv_]]:= iv;
	


End[ (* PrivateDiffAccess *) ];


Clear[DiffusionPrint, diffPrint];

diffPrint[x_,w_,f_,g_,i_] := diffPrint[x,w,f,g,i,x]

diffPrint[x_,w_,f_,g_,i_,alias_] := 
	Module[{xt = SequenceForm[alias,Subscript["t"]],
		    x0 = SequenceForm[alias,Subscript["0"]],
		    ft,gt},
		ft = f /. x -> xt;
		gt = g /. x -> xt;
		If [ NumberQ[gt] && gt == 0,
			Print["d",xt," = ",ft," dt","  ; ",x0," = ",i],
		(* Else non trivial dispersion *)
			Print["d",xt," = ",
				  ft," dt + ",gt," d",w,Subscript["t"],"  ; ",
				  x0," = ",i]
		];
	];


DiffusionPrint[diffusion[lab_,lab_,f_,g_,i_]] := Module \
	[	{},
		Print[lab, Subscript["t"],
			" = Weiner Process with scale ",g];
	]

DiffusionPrint[diffusion[x_,wp_,f_,g_,i_]] := 
	diffPrint[x,wp,f,g,i];
	
DiffusionPrint[diffusion[x_,wp_,f_,g_,i_],alias_String] :=
	diffPrint[x,wp,f,g,i,alias];


Clear[DiffusionSimulate]

Begin["PrivateDiffSimulate`"];


Clear[nRandom]
	nRandom[sigma2_] := 		(* single N[0,sigma2] *)
		N[ Sqrt[-sigma2 2 Log[ Random[] ]] Cos[Random[Real,{0,2 N[Pi]}]]];

	nRandom[sigma2_, n_] := 	(* list of n N[0,sigma2] *)
		N[ Sqrt[-sigma2 2 Log[Table[Random[],{n}]]] Cos[Table[Random[Real,{0, 2 N[Pi]}],{n}]]];
	
DiffusionSimulate::symbolic = 
	"This diffusion contains non-numeric parameters.";

DiffusionSimulate[diffusion[w_,w_,_,_,iv_], n_, dt_:0.01] := Module \
[	{i},
	Transpose[{ Table[i*dt,{i,0,n-1}],
			  FoldList[Plus,iv,nRandom[dt,n-1]] }
	]
]
	
DiffusionSimulate[d_diffusion, n_, dt_:0.01] := Module		\
[	{x,mu,sigma,sym,x0, mQ,sQ,xQ, dW, i=0},
	nextValue[{time_,current_}] := Module	\
	[	{dm,db},
		i  = i + 1;
		dm = mu[current,time] dt;
		db = sigma[current,time] dW[[i]];
		(* Print["dt ",dt," db ",db]; *)
		{time+dt,current+dm+db} // N
	];
	mu    = DiffusionDrift[d,Function];
	sigma = DiffusionDispersion[d,Function];
	x0    = N[DiffusionInitialValue[d]];
	mu11    = mu[1,1];	(* simple test for symbolic presence *)
	sigma11 = sigma[1,1];
	{mQ,sQ,xQ} = {NumberQ[mu11], NumberQ[sigma11], NumberQ[x0]};
	If[ mQ && sQ && xQ,
		 dW    = nRandom[dt,n];
		 NestList[nextValue,{0,x0},n-1],
	(* Else *)
		 Message[DiffusionSimulate::symbolic];
		 Print["mu[1,1] = ", mu11, ";  sigma[1,1] = ",sigma11,
		 	   ";  x0 = ",x0];
	]
];


End[ (* "PrivateDiffSimulate`" *) ];


Clear[A, DiffusionA, OperatorA];

DiffusionA[expr_] := A[expr];

A[expr_] := OperatorA[findDiffusion[expr]][expr];

OperatorA[d_diffusion][expr_] := Module	\
	[ {drift,disp, fd,sd, ggp,dt,xt,e=expr},
		e = expr /. d -> xt;
		drift = replaceDiffusion[DiffusionDrift[d],d,xt];
		disp  = replaceDiffusion[DiffusionDispersion[d],d,xt];
		fd  = D[e,xt];
		sd  = D[fd,xt];
		(* Print ["Drift ------>", drift," . ", fd]; *)
		(* Print ["Dispersion-->", disp," . ",disp,"' . ",sd]; *)
		drift= drift fd;
		disp = 1/2 disp disp sd;
		dt = D[e,Global`t];
		(* Print["Final three components"]; *)
		(* Print[drift, " + " , disp, " + ", dt]; *)
		Simplify[drift + disp + dt] /. xt -> DiffusionSymbol[d]
	];


Clear[Ito, OperatorIto, DiffusionIto, ItoInvert, ItoSymbol];
	
Options[Ito] = {ItoInvert->True, ItoSymbol->Z};

Ito[expr_, opts___Rule] := OperatorIto[findDiffusion[expr]][expr, opts];

DiffusionIto[expr_,opts___Rule] := Ito[expr,opts];

OperatorIto[d_diffusion][expr_, opts___Rule] := Module	\
	[ {drift, disp, e, xt, dSymbol, wp, invRules, ns, invert, start},
		dSymbol = DiffusionSymbol[d];
		wp = DiffusionWeinerProcess[d];
		e = replaceDiffusion[expr,d,xt];
		drift = OperatorA[d][expr] /. dSymbol -> xt;
		disp  = D[e,xt] (DiffusionDispersion[d] /. dSymbol -> xt);
		(* Print["Drift with xt ", drift,";   Disp with xt ", disp]; *)
		Print["d",e /. xt->dSymbol," = (",
				drift /. xt->dSymbol,") dt + (",
				disp /. xt->dSymbol,
				") d",DiffusionSymbol[wp],Subscript[t]];
		invert = ItoInvert /. {opts}  /. Options[Ito]; 
		If[invert,
			ns = ItoSymbol /. {opts} /. Options[Ito];
			(* Print["Solving ", ns, "  ==  ", e/.xt -> dSymbol]; *)
			invRules = First[Solve[ns==e,xt]];
			Print["Inversion rule... ", invRules /. xt -> dSymbol];
			(* Substitute in the drift and dispersion functions *)
			drift = drift /. invRules; drift=Simplify[drift];
			disp  = disp /. invRules; disp=Simplify[disp],
		(* Else *)
			ns = Head[expr];
		];
		(* Fix the starting location *)
		start = e/.xt->DiffusionInitialValue[d];
		(* Make the new diffusion in the global context *)
		If[ns===Z, Unprotect[Z]; Clear[Z]];
		DiffusionMake[ns,wp,drift/.xt->dSymbol,disp/.xt->dSymbol,start];
		Protect[Z];
		ns
	];


Clear[FeynmanKac, T, Infinitesimal];

DiffusionFeynmanKac::usage = "DiffusionFeynmanKac[pde, bv] is an
	alias for Feynman-Kac.";
	
FeynmanKac::usage="FeynmanKac[pde, bv] attempts to solve the pde
	represented by the first argument using the Feynman-Kac
	theorem.  The result, if successful, is an expectation over
	a diffusion.  Solution assumes interval [0,T].";
	
FeynmanKac::fails = "Could not solve this PDE via Feynman-Kac.";


Clear[InfRules]

InfRules = 	\
{	D[f_[x_,t_],t_]+m_.D[f_[x_,t_],x_]+s_.D[f_[x_,t_],{x_,2}]:>
			Infinitesimal[m,PowerExpand[Sqrt[2 s]]][ f[x,t] ],
	D[f_[x_,t_],t_] +                  s_.D[f_[x_,t_],{x_,2}]:>
			Infinitesimal[0,PowerExpand[Sqrt[2 s]]][ f[x,t] ],
	D[f_[x_,t_],t_] + m_. D[f_[x_,t_],x_]                    ->
			Infinitesimal[m,     0   ][ f[x,t] ]
};


Clear[FKW, FeynmanKac];			WeinerProcessMake[FKW];

FKRules :=		\
{
Infinitesimal[m_,s_][f_[x_,t_]] - rho_[x_,t_] f_[x_,t_] + u_[x_,t_] ->
				FK[f,diffusion[x,HoldForm[W],m,s,0],rho[x,t],u[x,t]],
Infinitesimal[m_,s_][f_[x_,t_]] - rho_. f_[x_,t_] + u_[x_,t_] ->
				FK[f,diffusion[x,HoldForm[W],m,s,0],rho,u[x,t]],
Infinitesimal[m_,s_][f_[x_,t_]] - rho_[x_,t_] f_[x_,t_] + u_. ->
				FK[f,diffusion[x,HoldForm[W],m,s,0],rho[x,t],u],
Infinitesimal[m_,s_][f_[x_,t_]] - rho_. f_[x_,t_] + u_. ->
				FK[f,diffusion[x,HoldForm[W],m,s,0],rho,u],
Infinitesimal[m_,s_][f_[x_,t_]] - rho_. f_[x_,t_] ->
				FK[f,diffusion[x,HoldForm[W],m,s,0],rho,0],
Infinitesimal[m_,s_][f_[x_,t_]] ->
				FK[f,diffusion[x,HoldForm[W],m,s,0],0,0]
};

DiffusionFeynmanKac[pde_, bv_] := FeynmanKac[pde,bv];

FeynmanKac[pde_,bv_] :=  Module		\
	[ {step1, step2, f,d,r,u},
		step1 = pde   /. InfRules;
		step2 = step1 /. FKRules;
		If [ Head[step2] === FK,
			{f,d,r,u} = List @@ step2;
			Print["f = ",f,";  rho = ",r,";  u = ",u];
			d[[5]] = x;
			DiffusionPrint[d,"X"];
			Ave[Exp[-r (Global`T-t)] bv[HoldForm[X][Global`T-t]] ],
		(* Else *)
			Message [FeynmanKac::fails];
			{pde,bv}
		]
	];
