
\section*{Appendix B: Simulation in Mathematica}

\begin{verbatim}
(*
  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),
           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)] }
    ]

\end{verbatim}
