Deutsches R-Forum

Verfasst: **Mi Mai 12, 2021 10:55 am**

Hallo

Ich möchte in R mir die kritischen Werte für den KPSS Test ausrechnen, hänge aber an der Berechnung für die Teststatistik fest. Vielleicht kann mir ja einer von euch helfen.

Hier erstmal die Schritte, die ich bisher gemacht habe.

Code: Alles auswählen

y <- c()
y[1] <- 0
r <- c()
r[1] <- 1
  
    Length <- 250  
for ( t in 2:Length){
  r[t] <- (r[t-1] + rnorm(1,0,0))}
for (t in 1 : Length){
  y[t] <- r[t]  + rnorm(1, 0,10)
}

Ich simuliere eine Zeitreihe ohne konstante, die stationär ist. Im nächsten Schritt brauch ich die Residuen der Regression, um sowohl den Zähler als auch den Nenner zu berechnen.

Code: Alles auswählen

reg <- lm(y ~ c(1:Length) -1)

Die Berechnung des Zählers habe ich so gemacht:

Code: Alles auswählen

sumres <- c()
for (i in 1:T){
  sumres[i] <- sum((reg$residuals)[1:i])}
for( t in 1:T){
  sqsumres <- sum(sumres[i]^2)

Habe ich das richtig gemacht, oder sieht hier jemand einen Denkfehler von mir?
Leider weiß ich weitergehend nicht, wie ich den HAC estimator der Varianz von den Residuen bestimme. Hierzu hätte ich NeweyWest() aus dem sandwich package genommen. Allerdings ist mir nicht schlüssig wie ich damit auf die Varianz der Residuen komme.

Ihr wärd mir eine Riesenhilfe, wenn mir einer Tipps geben kann oder mich auf Fehler aufmerksam macht.

Vielen Dank und beste Grüße,
Lukas

Verfasst: **Mi Mai 12, 2021 11:47 am**

Hallo Lukas,

leider kann ich Dir mit dem Kern Deines Problems nicht helfen. Trotzdem möchte ich Dich gerne auf Kleinigkeiten und auf die Funktion cumsum aufmerksam machen. Du hast programmiert:

Code: Alles auswählen

r <- c()
r[1] <- 1
  
    Length <- 250  
for ( t in 2:Length){
  r[t] <- (r[t-1] + rnorm(1,0,0))}

In der ersten Zeile teilst Du R mit, dass Du einen Vektor brauchen wirst, der r heißt. Alternativ kannst Du R auch mitteilen, dass Du einen Vektor brauchst, der numerische Werte enthalten und 250 Beobachtungen fassen soll. Das hat den Vorteil, dass R gleich festlegen kann, wieviel Speicher es dafür vom Betriebssystem anfordern muss und unnötige weitere Kommunikation mit dem Betriebssystem entfällt. Ist bei 250 egal, wird bei großen Zahlen aber ein wichtiger Weg, Rechnungen zu beschleunigen:

Code: Alles auswählen

r <- numeric(250)
r[1] <- 1

Alternativ dazu würde ich jetzt r gleich definieren, vorläufig als einen Vektor, der mit 1 anfängt, gefolgt von 250 Zufallszahlen:

Code: Alles auswählen

r <- c(1, rnorm(249))

Ich weiß, das ist jetzt inhaltlich noch falsch, denn wir wollen ja nicht 249 Zufallswerte sondern die kummulierten Summen dieser Zufallszahlen. Und genau das macht die Funktion cumsum. Probier mal:

Code: Alles auswählen

cumsum(c(1, 1, 1, 1, 1, 1, 1, 1, 1))
cumsum(c(1, -1, 1, -1, 1, -1, 1, -1))

Damit wird aus den eingangs zitierten 5 Zeilen Code jetzt

Code: Alles auswählen

r <- cumsum(c(1, rnorm(249)))

Das beschreibt nicht nur das Problem viel kompakter sondern wird auch wesentlich schneller (wenn man dann irgendwann nicht mehr 250 sondern vieleicht 250 Mio Durchgänge rechnen will).

Zur Kontrolle, ob diese Schreibweise wirklich random walks produziert kannst Du Dir gleich mehrere davon plotten lassen:

Code: Alles auswählen

par(mfrow=c(3,3), mar = c(2,2,1,1))
replicate(9, {r <- cumsum(c(1, rnorm(249)));plot(r, type = "l")})

HTH,
Bernhard

Verfasst: **Mi Mai 12, 2021 11:55 am**

Hi Bernhard,
vielen Dank.
Ich merke oft, dass ich in Sachen effizienz noch einiges verbessern muss, von daher bin ich für solche Tipps auch immer sehr dankbar

Ich hoffe, dass ich bald meine anderen Probleme gelöst habe, sodass ich mich dann der Verbesserung meines Codes widmem kann.

Danke und beste Grüße,
Lukas

Verfasst: **Mi Mai 12, 2021 1:28 pm**

Gude!

Was ist mit dem KPSS-Test aus der library tseries?
https://www.rdocumentation.org/packages ... /kpss.test

Dank&Gruß
Schubbiaschwilli

Verfasst: **Mi Mai 12, 2021 1:35 pm**

Servus,
ich weiß, dass man das da bekommt. Ich wollte das nur mal selber nachrechnen

.

Danke und beste Grüße,
Lukas

Verfasst: **Mi Mai 12, 2021 1:43 pm**

Das aber natürlich trotzdem praktisch, wenn man sowas fertig findet, dann kann man ggf. in deren open-source Code nachschauen, wie die das Problem gelöst haben. Siehe hier: https://github.com/cran/tseries/blob/master/R/test.R ab Zeile 720.

LG,
Bernhard

Verfasst: **Mi Mai 12, 2021 1:47 pm**

Hi,
ja, genau. Vor allem um zu überprüfen, ob ichs richtig mache.
Danke für den Link, daran hab ich noch gar nicht gedacht

Danke und beste Grüße,
Lukas

Verfasst: **Mi Mai 12, 2021 2:06 pm**

Gude!

Man kann auch einfach den Namen der Funktion eingeben - Oder wäre das zu einfach?
Und sich dann eine MyOwn.kpss.test-Funktion bauen.

Dank&Gruß
Schubbiaschwilli

Code: Alles auswählen

> kpss.test
function (x, null = c("Level", "Trend"), lshort = TRUE) 
{
    if ((NCOL(x) > 1) || is.data.frame(x)) 
        stop("x is not a vector or univariate time series")
    DNAME <- deparse(substitute(x))
    null <- match.arg(null)
    x <- as.vector(x, mode = "double")
    n <- length(x)
    if (null == "Trend") {
        t <- 1:n
        e <- residuals(lm(x ~ t))
        table <- c(0.216, 0.176, 0.146, 0.119)
    }
    else if (null == "Level") {
        e <- residuals(lm(x ~ 1))
        table <- c(0.739, 0.574, 0.463, 0.347)
    }
    tablep <- c(0.01, 0.025, 0.05, 0.1)
    s <- cumsum(e)
    eta <- sum(s^2)/(n^2)
    s2 <- sum(e^2)/n
    if (lshort) 
        l <- trunc(4 * (n/100)^0.25)
    else l <- trunc(12 * (n/100)^0.25)
    s2 <- .C(tseries_pp_sum, as.vector(e, mode = "double"), 
        as.integer(n), as.integer(l), s2 = as.double(s2))$s2
    STAT <- eta/s2
    PVAL <- approx(table, tablep, STAT, rule = 2)$y
    if (!is.na(STAT) && is.na(approx(table, tablep, STAT, rule = 1)$y)) 
        if (PVAL == min(tablep)) 
            warning("p-value smaller than printed p-value")
        else warning("p-value greater than printed p-value")
    PARAMETER <- l
    METHOD <- paste("KPSS Test for", null, "Stationarity")
    names(STAT) <- paste("KPSS", null)
    names(PARAMETER) <- "Truncation lag parameter"
    structure(list(statistic = STAT, parameter = PARAMETER, p.value = PVAL, 
        method = METHOD, data.name = DNAME), class = "htest")
}
<bytecode: 0x00000151ecc2af60>
<environment: namespace:tseries>
>

Verfasst: **Mi Mai 12, 2021 2:12 pm**

hi,
ja ich hab noch nicht daran gedacht mir den Code für die funktion von kpss.test anzugucken.
Werde ich mir merken für die Zukunft

Danke und beste Grüße,
Lukas

Verfasst: **Mi Mai 12, 2021 3:27 pm**

schubbiaschwilli hat geschrieben: ↑Mi Mai 12, 2021 2:06 pmMan kann auch einfach den Namen der Funktion eingeben - Oder wäre das zu einfach?

Sieht bei mir so aus:

Code: Alles auswählen

> library(tseries)
Fehler in library(tseries) : es gibt kein Paket namens ‘tseries’

deshalb der Blick auf GitHub. Ich lad mir nicht mehr jedes Package herunter, bloß weil sich eine Frage im Forum darauf bezieht.

LG,
Bernhard

Deutsches R-Forum

KPSS - Test

KPSS - Test

Re: KPSS - Test

Re: KPSS - Test

Re: KPSS - Test

Re: KPSS - Test

Re: KPSS - Test

Re: KPSS - Test

Re: KPSS - Test

Re: KPSS - Test

Re: KPSS - Test