Monte-Carlo-Simulation der Korrelation zwischen zwei Brownsche Bewegung (continuous random walk)

Question

y <- cumsum(rnorm(100,0,1)) # random normal, with small (1.0) drift. 
y.ts <- ts(y) 
x <- cumsum(rnorm(100,0,1)) 
x 
x.ts <- ts(x) 
ts.plot(y.ts,ty= "l", x.ts) # plot the two random walks 


Regression.Q1 = lm(y~x) ; summary(lm2) 
summary(Regression.Q1) 

t.test1 <- (summary(Regression.Q1)$coef[2,3]) # T-test computation 


y[t] = y[t-1] + epsilon[t] 
epsilon[t] ~ N(0,1) 
set.seed(1) 
t=1000 
epsilon=sample(c(-1,1), t, replace = 1) # Generate k random walks across time {0, 1, ... , T} 


N=T=1e3 
y=t(apply(matrix(sample(c(-1,1),N*T,rep=TRUE),ncol=T),1,cumsum)) 
y[1]<-0 
for (i in 2:t) { 
    y[i]<-y[i-1]+epsilon[i] 
} 

Ich muß:

den Vorgang wiederholen 1.000mal (Monte-Carlo-Simulationen), nämlich Aufbau eine Schleife um das vorherige Programm und jedes Mal, Speichere die t Statistiken. Sie haben eine Sequenz von 1 000 t-Tests: S = (t-test1, t-test2, ..., t-test1000). Zählen Sie die Anzahl der Zeit, die der absolute Wert der 1.000 t-Tests> 1,96 ist, der kritische Wert bei einem Signifikanzniveau von 5%. Wenn die Serie I (0) wäre, hätten Sie ungefähr 5% gefunden. Dies wird hier nicht der Fall sein (falsche Regression).

Was muss ich hinzufügen, um die entsprechenden Koeffizienten zu speichern?

Answer 1

Ihr geposteter Code bezogen auf y[t] = y[t-1] + epsilon[t] ist nicht wirklich funktionierender Code, aber ich kann sehen, dass Sie versuchen, alle 1000 * 2 Random Walk zu speichern. Eigentlich ist das nicht nötig. Wir kümmern uns nur um den T-Score und nicht darum, was diese Realisierungen von Random Walk sind.

Für diese Art von Problem, bei dem wir versuchen, eine Prozedur viele Male zu replizieren, ist es praktisch, zuerst eine Funktion zu schreiben, um eine solche Prozedur einmal auszuführen. Sie hatten bereits einen guten Arbeitscode dafür; wir müssen es nur in einer Funktion wickeln (diesen unnötigen Teil wie plot entfernen):

sim <- function() { 
    y <- cumsum(rnorm(100,0,1)) 
    x <- cumsum(rnorm(100,0,1)) 
    coef(summary(lm(y ~ x)))[2,3] 
    } 

Diese Funktion keine Eingabe erfolgt; Es gibt nur den T-Wert für ein Experiment zurück.

Jetzt werden wir das 1000 Mal wiederholen. Wir können eine for Schleife schreiben, aber Funktion replicate ist einfacher (lesen ?replicate wenn nötig)

S <- replicate(1000, sim()) 

Hinweis: Das wird einige Zeit dauern, viel langsamer, als es für eine so einfache Aufgabe sein sollte, weil beide lm und summary.lm sind langsam. Ein viel schnellerer Weg wird später gezeigt.

Jetzt ist S Vektor mit 1000 Werten, das ist die "eine Folge von 1000 t-Tests" Sie wollen. Um „die Anzahl der Zeit der absolute Wert der 1000 t-Tests> 1,96“, können wir verwenden nur

sum(abs(S) > 1.96) 
# [1] 756 

Das Ergebnis 756 ist genau das, was ich bekommen; Sie werden etwas anderes bekommen, da die Simulation zufällig ist. Aber es wird immer eine ziemlich große Anzahl sein, wie erwartet.

---

Eine schnellere Version von sim:

fast_sim <- function() { 
    y <- cumsum(rnorm(100,0,1)) 
    x <- cumsum(rnorm(100,0,1)) 
    y <- y - mean(y) 
    x <- x - mean(x) 
    xty <- crossprod(x,y)[1] 
    xtx <- crossprod(x)[1] 
    b <- xty/xtx 
    sigma <- sqrt(sum((y - x * b)^2)/98) 
    b * sqrt(xtx) * sigma 
    } 

Diese Funktion berechnet einfache lineare Regression ohne lm, und T-Score ohne summary.lm.

S <- replicate(1000, fast_sim()) 
sum(abs(S) > 1.96) 
# [1] 778 

Eine alternative Möglichkeit ist cor.test zu verwenden:

fast_sim2 <- function() { 
    y <- cumsum(rnorm(100,0,1)) 
    x <- cumsum(rnorm(100,0,1)) 
    unname(cor.test(x, y)[[1]]) 
    } 

S <- replicate(1000, fast_sim()) 
sum(abs(S) > 1.96) 
# [1] 775 

ist eine Benchmark Let haben:

system.time(replicate(1000, sim())) 
# user system elapsed 
# 1.860 0.004 1.867 

system.time(replicate(1000, fast_sim())) 
# user system elapsed 
# 0.088 0.000 0.090 

system.time(replicate(1000, fast_sim2())) 
# user system elapsed 
# 0.308 0.004 0.312 

cor.test ist viel schneller als lm + summary.lm, aber manuelle Berechnung ist noch schneller!