(*
   Burton Rosenberg
   Estimating Model Parameters in Markets with Collars
   11/27/2000

   Added Black-Scholes equations

*)

(* include packages *)
<< Statistics`ContinuousDistributions`
<< Statistics`DescriptiveStatistics`

(* trace generation routines *)

weinerWalk[ n_, m_, sd_ ] :=
    Module [
        { s = 0,
            sl = {0}
        },
        Do [
            s = s + Random[NormalDistribution[m, sd]] ;
            sl = Append[sl, s]
            , {n}
        ] ;
        sl 
    ]
    
collarX[x_, s_, l_, u_ ] := 
    If[ ( x > (s+u) ), (s+u) , 
       If[ (x < (s-l) ),  (s-l), x
       ]
    ]

applyCollar[ walk_, l_, u_ ] :=
    Module [
        { walkL = {},
          x = walk[[1]]
        }, 
        Do [
            x = collarX[ walk[[i]], x, l, u ] ;
            walkL = Append[ walkL, x ]
            , {i, 1, Length[walk]}  
        ] ;
        walkL
    ]

classifyDays[ a_, b_ ] :=
    Module [
        { c = {} },
        Do [ 
            If [ ( a[[i]]==b[[i]] ),
                c = Append[ c, 0 ],
                c = Append[ c, 1 ]
            ]
            , {i, 1, Length[a]}
        ] ;
        c
    ]

generateTrace[ n_, mean_, sd_, lower_, upper_] :=
    Module [
        { rw, cw, cd },
        rw =  weinerWalk[ n, mean, sd ] ;
        cw = applyCollar[ rw, lower, upper ] ;
        cd = classifyDays[ rw, cw ] ;
        (* return [ ideal walk, observed walk, days classified } *)
        (* codes for days classified: 0==at-market, 1==at-limit *)
        { rw, cw, cd }
    ]


(* end trace generation routines *)

(* helper functions to estimation routines *)

countTransitions[ daysClassified_ ] :=
    Module [
        { result = { 0, 0, 0, 0 }, 
          index
        },
        Do [
          If[ ( daysClassified[[i]]==0 ), 
              index = 1, index = 3 ] ;
          If[ ( daysClassified[[i+1]]==0 ),
              index = index, index = index+1 ] ;
          result[[index]]++
          , {i, 1, Length[daysClassified]-1}
        ] ;
        (* { marketMarket, marketLimit,
             limitLimit, limitMarket } *)
        result
    ]

waitingTimes[ daysClassified_ ] :=
    Module [
       { waitT = {},
         classOfPreviousDay,
         currentWait = 0
       },
       (* assert: daysClassified[[1]]==0 *)
       Do [
           If[ (daysClassified[[i]]==0),
               currentWait++ ;
               , (* else *)
               If [ (classOfPreviousDay==0),
                   waitT = Append[ waitT, currentWait ] ;
                   currentWait = 0
               ]
           ] ;
           classOfPreviousDay = daysClassified[[i]] 
           , {i, 1, Length[daysClassified]}
       ] ;
       waitT
    ]

lastMarketDay[ daysClassified_ ] :=
    Module [ { lmd },
        (* assert: daysClassified[[1]]==0 *)
        Do [
            If [ (daysClassified[[i]]==0),
                lmd = i
            ]
            , {i, 1, Length[daysClassified] }
        ] ;
        lmd
    ]


f[ mu_, sigma_, ll_, ul_ ] := 
    1 - Integrate[
           ( 1 / (sigma Sqrt[ 2 Pi ]) ) Exp[ - (t-mu)^2 / ( 2 sigma^2 ) ], 
           { t, ll, ul }
        ]

(* estimation routines *)

estimateMean[ aTrace_ ] :=
    Module[ { lmd = lastMarketDay[ aTrace[[3]] ] },
       (*note: arrays count from 1, differs from text *)
       (aTrace[[1]])[[lmd]]/(lmd-1)
    ]

waitingTimeMethod[ aTrace_, lowerCollar_, upperCollar_ ] :=
    Module[ { wt, p },
       wt = waitingTimes[ aTrace[[3]] ] ;
      (* if no limit days, create one waiting time, with value beyound maximum *)
       If [ (Length[wt]==0), 
(*Print["waitingTimeMethod: no limit days - append maximum possible wait time"] ;*)
           wt = { Length[aTrace[[1]]] } ] ; 
       p = 1 / Mean[ wt ] ;
       x /. FindRoot[ f[estimateMean[ aTrace ], x,
           -lowerCollar, upperCollar]==p, {x,1},
           MaxIterations->100, AccuracyGoal->15 ]//N
    ]

transitionCountingMethod[ aTrace_, lowerCollar_, upperCollar_ ] :=
    Module[ { cx, p },
        cx = countTransitions[ aTrace[[3]] ] ;
        p = cx[[2]]/(cx[[1]]+cx[[2]]) ;
If[ (p==0), Print["transitionCountingMethod: no limit days"] ] ;
        x /. FindRoot[ f[estimateMean[ aTrace ], x, 
            -lowerCollar, upperCollar]==p, {x,1},
            MaxIterations->100, AccuracyGoal->15 ]//N
    ]

generalizedMomentsMethod[ aTrace_ ] :=
    Module [ 
        {   data = aTrace[[1]],
            dayType = aTrace[[3]],
            mu = estimateMean[ aTrace ],
            listLength = Length[ aTrace[[1]] ],
            sum, i, j 
        } ,
        (* loop invariant, sum contains all data up to day i;
           day i is a market day *)
        sum = 0 ;
        i = 1 ;
        While [ (True),
           j = i+1 ;
           While [ (j<=listLength && dayType[[j]]==1), j++ ] ;
           If [ (j>listLength), Break[] ] ;
           sum += ((data[[j]]-data[[i]])-(j-i)*mu)^2 ;
           i = j ;
        ] ;
        Sqrt[sum/listLength]
    ]

observedMoments[ aTrace_ ] :=
    Module [
        {  data = aTrace[[2]],
           mu, sum
        },
        sum = 0 ;
        mu = data[[Length[data]]]/(Length[data]-1) ;
        Do [
            sum += ( data[[i+1]]-data[[i]]-mu )^2 ;
            ,{ i, 1, Length[data]-1 }
        ] ;
        Sqrt[sum/(Length[data]-2)]
    ]

(* Black Scholes *)

blackScholesDs[ strikePrice_, marketPrice_, timeToExpiration_,
    volatility_, riskFreeRate_ ] :=
    Module [
        {
            a = Log[marketPrice/strikePrice],
            b1 = timeToExpiration ( riskFreeRate + volatility^2/2 ),
            b2 = timeToExpiration ( riskFreeRate - volatility^2/2 ),
            c = volatility Sqrt[timeToExpiration]
        },
        { (a+b1)/c, (a+b2)/c }
    ]
   

(* volatility of 48% is input as .48. Interest rate of 6% is input as .06.
   n days is input as n/360 *)

blackScholes[ strikePrice_, marketPrice_, timeToExpiration_,
    volatility_ , riskFreeRate_] :=
    Module [
        { d = blackScholesDs[ strikePrice, marketPrice, timeToExpiration,
              volatility, riskFreeRate ] 
        },
        marketPrice CDF[NormalDistribution[0,1],d[[1]]] -
            strikePrice Exp[- riskFreeRate timeToExpiration] CDF[NormalDistribution[0,1],d[[2]]]
    ]

(* compare methods *)

pointIt[ x_,ys_ ] :=
    Table[ { x, ys[[n]] }, { n,Length[ys] } ]

boxIt[ x_, ys_ ] :=
    Module [ { q1 = Quantile[ys,1/4],
               q2 = Quantile[ys,1/2],
               q3 = Quantile[ys,3/4],
               d = .1
              },
              Line[ {{x-d,q2},{x-d,q3},{x+d,q3},{x+d,q1},{x-d,q1},
                     {x-d,q2},{x+d,q2}} ] 
    ]

compete[ trials_, collar_ ] :=
    Module [ { tr, n = 250, m=0, sd=1, lc=collar, uc=collar,
         wtr = {}, gmr = {}, om = {} },
       Do [
           tr = generateTrace[ n, m, sd, lc, uc ] ;
           wtr = Append[ wtr, waitingTimeMethod[ tr, lc, uc ] ] ;
           gmr = Append[ gmr, generalizedMomentsMethod[ tr ] ] ;
           om = Append[ om, observedMoments[ tr ] ]
           , { i, 1, trials }
       ] ;
       ListPlot[ Join[ pointIt[1, wtr], pointIt[2, gmr], pointIt[3, om] ],
          PlotRange->{{0,4}, Automatic}, Ticks->{None,Automatic},
          Epilog->{ RGBColor[1,0,0], boxIt[1,wtr], boxIt[2,gmr], boxIt[3,om] } 
       ] 
       (* ;{ Quantile[wtr,1/2], Quantile[gmr,1/2] } *)
    ]

(* price options *)

(* Antithetic walks, for faster convergence in simulation.
   See p 89, Implementing Derivatives Modesl, Les Clewlow and
   Chris Strickland, Wiley 1998. HG6024.A3C584 (ISBN 0-471-96651-7)
   *)
   
antitheticWalks[ n_, m_, sd_, lower_, upper_ ] :=
    Module [
        { sp=0, sn=0, spl={}, snl={} },
        Do [
            sp += m + sd * Random[NormalDistribution[0,1]] ;
            spl = Append[ spl, sp ] ;
            sn += m - sd * Random[NormalDistribution[0,1]] ;
            snl = Append[ snl, sn ] ;
            , {n}
        ] ;
        { spl, snl, applyCollar[ spl, lower, upper ], applyCollar[ snl, lower, upper ] }
    ]

priceOptions[ trials_, collar_ ] :=
    Module [
        {
           tr, winnings1=0, winnings2=0,
           lc, uc, n=30, sd=1/30 (* yearly volatility == (sd=1/30)/Sqrt[1/360] = .632456 *),
           m, 
           rfr = 0
        },
        m = rfr - sd^2/2 ;
        Do [
            lc = collar * sd ;
            uc = lc ;
            tr = antitheticWalks[ n, m, sd, lc, uc ] ;
            If[ ( (tr[[1]])[[n-1]]>0 ),
                  winnings1 += (Exp[(tr[[1]])[[n-1]]]-1)/2
            ] ;
            If[ ( (tr[[2]])[[n-1]]>0 ),
                  winnings1 += (Exp[(tr[[2]])[[n-1]]]-1)/2
            ] ;
            If[ ( (tr[[3]])[[n-1]]>0 ),
                  winnings2 += Exp[(tr[[3]])[[n-1]]]-1
            ]
            , { i, 1, trials }
        ] ;
        { Exp[- n rfr /360] winnings1/trials, Exp[- n rfr /360 ] winnings2/trials }
    ]

(* Draw price charts *)

ro[l_] := Module[ { ltemp = {} }, 
                  Do[ ltemp = Append[ltemp, {i, Exp[l[[i]]]}], 
                      {i, 1, Length[l]}]; 
                  ltemp
                ]

ro2[l1_,l2_, d_, s_] :=
    Module [
        { lt={} },
        Do [ lt = Append[ lt, 
                   {i, s (Exp[l1[[i]]]-Exp[l2[[i]]])-d}
                   ] 
            , {i, 1, Length[l1]}
            ] ;
        lt
    ]

myShow2[l1_,l2_, d_, s_] := Show[
                       Graphics[{RGBColor[1, 0, 0], 
                                Line[ ro[l1] ], 
                                RGBColor[0, 0, 1], 
                                Line[ ro[l2] ],
                                RGBColor[0,0,0],
                                Line[ ro2[l1,l2,d,s]]}
                       ],
                       {Axes -> Automatic}]
                       
myShow[l1_,l2_, d_, s_] := Show[
                       Graphics[{RGBColor[1, 0, 0], 
                                Line[ ro[l1] ], 
                                RGBColor[0, 0, 1], 
                                Line[ ro[l2] ]}
                       ],
                       {Axes -> Automatic}]

priceCharts[ c_ ] :=
    Module [
        { tr, m=0, sd=1/40, n=50, lc, uc  },
        lc = sd * c ;
        uc = sd * c ;
        tr = generateTrace[ n, m, sd, lc, uc ] ;
        myShow[ tr[[1]], tr[[2]], -1/2, 2 ]
    ]
    

(* ***NOTES*** *)

(* after running compete, use Display["filename.ps", %, "EPS"] to save image as PS *)
