Simulation eines stochastischen Integrals

7

Ich versuche, Übung 3.9.10 auf S. 22 zu lösen . 66 von Ubbo F. Wiersemas "Brownian Motion Calculus" (John Wiley & Sons, 2008) , in dem das stochastische Integral simuliert werden soll

01B.(t) dB.(t)
indem zunächst eine Partition von in Teilintervalle verwendet wird und Simulationen des diskreten stochastischen Integrals für dieses n , wobei dieser Vorgang wiederholt wird, indem n wiederholt auf 2 ^ {11} verdoppelt wird . In der Übung werden die Ergebnisse mit dem Mittelwert und der Varianz der geschlossenen Form des Integrals verglichen, nämlich \ int_0 ^ 1 B (t) \ dB (t) = \ frac {1} {2} B (1) ^ 2 - \ frac {1} {2}, die 0 bzw. \ frac {1} {2} sind (da B (1) eine normale Standardzufallsvariable ist).[0,1]]n=281000
ich(n)=ich=0n- -1B.(tich)(B.(tich+1)- -B.(tich))
nn211
01B.(t) dB.(t)=12B.(1)2- -12
012B.(1)

Ich habe eine Simulation in R geschrieben, wie unten gezeigt. Die aus dieser Simulation resultierenden Abweichungen sind jedoch

> v_sim
        8         9        10        11 
0.4771895 0.4304475 0.5260542 0.4664552

die nicht auf den erwarteten Wert zu konvergieren scheinen . Das gleiche Phänomen wird beobachtet, wenn der Keim auf und geändert wird und wenn die Anzahl der Unterintervalle auf auf erhöht wird . Was mache ich falsch?0,5232122fünfzehn

NSIMS <- 1000L # number of simulations of the integral for every "level"
LEVELS <- 8L:11L # level n corresponds to 2^n sample points between 0 and 1
SEED <- 1

set.seed(SEED)

sims <- data.frame(simid = 1L:NSIMS) # sims contains one colum per every level.
                                        # The column simid is a dummy column
                                        # intended to avoid an error message
                                        # when columns are added to sims
                                        # on the fly inside the following loop.
                                        # This column is deleted after the loop
                                        # ends.
for (n in LEVELS)
{
        nticks <- 2^n # nticks is the number of sample points between 0 and 1
        delta <- 1/nticks
        ticks <- seq(from = 0, to = 1, by = delta)
        std <- sqrt(delta)
        sim <- vector(mode = "numeric", length = NSIMS)
        for (j in 1L:NSIMS)
        {
                b <- cumsum(c(0, rnorm(nticks, sd = std))) # b is a simulated
                                                        # Brownian motion path
                                                        # sampled at the tick
                                                        # marks.
                integral <- 0
                for (i in 1L:(length(b) - 1L))
                {
                        integral <- integral + b[i]*(b[i + 1] - b[i])
                }
                sim[j] <- integral
        }
        sims[, as.character(n)] <- sim
}
sims$simid <- NULL

m_sim <- sapply(sims, mean)
v_sim <- sapply(sims, var)
Evan Aad
quelle
Haben Sie die Varianz bei jedem n berechnet? Sie werden nie genau 0,5 erhalten. Sie müssen also feststellen, ob Sie Ihre Nullhypothese von 0,5 ablehnen können.
Greg Petersen
@ GregPetersen: Ich erwarte nicht genau zu sehen 0,5, aber ich erwarte eine konsequente Konvergenz zu 0,5, z.B 0.4771895,0.48044750.51605420.4964552.
Evan Aad
2
Das Programm macht das, was es sollte, obwohl es sowohl einfacher als auch effizienter gestaltet werden kann: Verwenden Sie ein cumsumanstelle der iSchleife. Matrix / Apply könnte hier auch dem Datenrahmen / Sapply usw. vorgezogen werden. Dies ist eine sehr gute Übung. Tipps: Geben Sie klar an, was Sie erwarten und warum, finden Sie die gemeinsame Verteilung dern summierte Variablen sagen X.ich.
Yves

Antworten:

4

Sie haben zwei Abweichungsquellen von 0,5 und nehmen nur Änderungen an einer davon vor. Infolgedessen erwarten Sie nicht die von Ihnen vorgeschlagene Konvergenz. [Stattdessen würde man im kleinsten Fall eine leichte Verbesserung erwartenn Werte, die dann mit zunehmendem Wert im Rauschen verschwinden n des Weiteren.]

Die erste Quelle der Divergenz wird in Ihrer Frage erkannt - die Approximation des kontinuierlichen Prozesses durch eine diskrete Approximation.

Die zweite Quelle ist die zufällige Variation aufgrund der Simulation. Bei einer endlichen Anzahl von Simulationen (b) erwarten Sie, dass sich die Berechnung der Stichprobenvariation von dem zugrunde liegenden "wahren" Wert unterscheidet, den Sie mit einer bestimmten diskreten Näherung als erhalten würden b.

Um zu sehen, ob dies tatsächlich das Problem ist, können Sie eine Folge von versuchen n Werte, sagen wir n = 100.1000.10000 (und vielleicht mehr, wenn Sie können) bei jedem n. Alternativ können Sie eine theoretische Varianz der von Ihnen ausgedruckten Mengen (und damit einen Standardfehler) berechnen, anhand derer Sie feststellen können, ob Ihre Ergebnisse eindeutig nicht mit 1/2 übereinstimmen (asn wird groß Sie sollten feststellen, dass die Inkonsistenz schwerer zu erkennen ist, da sie viel größer erfordern würde b zu identifizieren).

Glen_b -State Monica
quelle
4

Lassen Xi:=B(ti)[B(ti+1)B.(tich)]] zum ich=0 zu n- -1. Die rv.s. X.ich sind gaußsch und unabhängig aufgrund der Unabhängigkeit der Inkremente von B.(t). AußerdemE.(X.ich)=0 und

Var(X.ich)=E.[X.ich2]]=E.[B.(tich)2]]E.[(B.(tich+1)- -B.(tich))2]]=tich(tich+1- -tich),
damit Var(X.ich)=ich/.n2 wenn die Punkte tich: =ich/.nwerden verwendet. Dann
Var(ich(n))=ich=0n- -1Var(X.ich)=1n2ich=0n- -1ich=1n2(n- -1)n2=[1- -1n]]12.
Bis zur Laufzeit 1/.n, Die Verteilung von ich(n) bleibt für alle fast gleich n. Sagen Sie, indem Sie eine Zahl zeichnenN.simvon unabhängig ich(n)erhalten wir eine Stichprobe der Größe N.sim von einer Verteilung abhängig nur geringfügig von n und das ist fast der normale mit Mittelwert 0 und Varianz 1/.2. Im Programm sind die v_simentsprechenden Elemente n(sofern sie nicht nklein sind) einfach die Varianzen unabhängiger StichprobenN.sim von Norm(0,1/.2).

Die Nummer nDie Anzahl der Punkte ist hier aufgrund des betrachteten spezifischen Integranden nahezu irrelevant. Die Situation kann mit der Bewertung des nichtstochastischen Integrals verglichen werden01tdtunter Verwendung der Trapezregel. Das Ergebnis ändert sich nicht, wenn mehr Trapeze verwendet werden. Offensichtlich wären die Dinge anders, wenn ein anderer angepasster Prozess oder eine andere angepasste Funktion als Integrand verwendet würde, ein dichteres Design[tich]]ich was dann zu einer kleineren Vorspannung führt.

Yves
quelle
2
Kleiner Tippfehler denke ich in der ersten Anzeigegleichung mit fehlenden Quadraten.
JDODS
0

Um die Konvergenz zu simulieren, simulieren Sie den Mittelwert des quadratischen UNTERSCHIEDS zwischen dem Integral basierend auf n Schritten und dem Integral basierend auf 2n Schritten. Dann das gleiche für die Differenz zwischen den Integralen für 2n Schritte und 4n Schritte. Dann der Unterschied zwischen 4n Schritten und 8n Schritten. Verwenden Sie für alle die gleiche Brownsche Bewegungsbahn (in der Zeiteinheit). Das gleiche Konvergenzmuster kann unter Verwendung der Varianz des UNTERSCHIEDS zwischen den oben erwähnten Integralpaaren gesehen werden. Abschnitt 3.4 meines Buches behandelt dies. Ubbo

user186735
quelle
-2

3.9.10 BdB.R.

Integral B.dB

Übung 3.9.10 wie folgt umformulieren:

Löschen Sie Zeile 5 und weiter bis zum Satz von Zeile 8, der mit Wiederholen beginnt

Löschen Sie auch den Satz, der mit Vergleichen beginnt

set.seed (123) # optional; positive Ganzzahlsimulationen <- 1000 # Anzahl der Simulationen

BdB.sim <- Funktion (n) {B.Schritte <- 2 * n # Anzahl der Zeitschritte im Brownschen Bewegungspfad

# integrals B.dB based on successively halved time partitions
I.n <-  numeric(sims)  # n steps
I.2n <- numeric(sims)  # 2*n steps

dt <- 1/B.steps
sqrt.dt <- sqrt(dt)
k.max <- 1+B.steps     # last timepoint index of B

# timepoints
# ---------------------------------------
k.seq.n <- seq(from=3,to=k.max,by=2)
length(k.seq.n)

k.seq.2n <- seq(from=2,to=k.max,by=1)
length(k.seq.2n)

k.seq.step <- seq(from=3,to=k.max-1,by=2)
length(k.seq.step)
# ---------------------------------------

for(sim.nr in 1:sims) {
    # each simulation starts by generating a Brownian Motion path of n timesteps, denoted B
    B <- c(0,cumsum(rnorm(B.steps,mean=0,sd=sqrt.dt)))  # Brownian Motion path timestep dt; B[1]=0
    for(k in k.seq.n) {I.n[sim.nr] <- I.n[sim.nr] + B[k-2]*(B[k]-B[k-2])}
    for(k in k.seq.2n) {I.2n[sim.nr] <- I.2n[sim.nr]  + B[k-1]*(B[k]-B[k-1])}
}

mean.sq.diff <- mean((I.n-I.2n)^2)
exact <- 1/(4*n)

return(data.frame(n=n,double=2*n,mean.sq.diff=mean.sq.diff,exact=exact))

}}

(n.seq <- c (2 ^ 8, 2 ^ 9, 2 ^ 10,2 ^ 11))

trail <- data.frame (n = NA, double = NA, mean.sq.diff = NA, genau = NA)

für (m in 1: Länge (n.seq)) {out <- BdB.sim (n.seq [m]) trail <- rbind (trail, out)}

Spur [-1,] # zeigt Konvergenz

write.csv (trail [-1,], "3.9.10 BdB convergence.csv", row.names = FALSE) # optional

--------------------------

Ubbo Wiersema 26. Februar 2018

user186735
quelle
Für aufeinanderfolgende Potenzen von n ist das Verhältnis 1 / (4 (n + 1)) über 1 / 4n gleich n / (n + 1), und da n eine Potenz von 2 ist, ist dieses Verhältnis gleich ½. Die mittlere quadratische Differenz verringert sich um die Hälfte, wenn die Anzahl der Zeitintervalle von 2ⁿ auf 2ⁿ⁺¹ erhöht wird.
user186735