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

*)

(* 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 }
        ] ;
        { mu, Sqrt[sum/(Length[data]-2)] }
    ]


(* sample runs *)

runAll[ n_, m_, sd_, lc_, uc_ ] :=
    Module[ { theTrace },
        theTrace = generateTrace[ n, m, sd, lc, uc ] ;
        { estimateMean[ theTrace ], 
          waitingTimeMethod[ theTrace, lc, uc ],  
          transitionCountingMethod[ theTrace, lc, uc ],
          generalizedMomentsMethod[ theTrace ],
          observedMoments[ theTrace ]
        }
    ]

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 = {} },
       Do [
           tr = generateTrace[ n, m, sd, lc, uc ] ;
           wtr = Append[ wtr, waitingTimeMethod[ tr, lc, uc ] ] ;
           gmr = Append[ gmr, generalizedMomentsMethod[ tr ] ]
           , { i, 1, trials }
       ] ;
       ListPlot[ Join[ pointIt[1, wtr], pointIt[2, gmr] ],
          PlotRange->{{0,3}, Automatic}, Ticks->{None,Automatic},
          Epilog->{ RGBColor[1,0,0], boxIt[1,wtr], boxIt[2,gmr] } 
       ] 
       (* ;{ Quantile[wtr,1/2], Quantile[gmr,1/2] } *)
    ]

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

(* \special(eps:c:\filename.ps x=1in y=3in} *)

