KPSS - Test

Methoden der Zeitreihenanalyse

Moderator: schubbiaschwilli

Antworten
Luk12Luk12
Beiträge: 25
Registriert: Di Dez 01, 2020 12:31 pm

KPSS - Test

Beitrag von Luk12Luk12 »

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.
Screenshot 2021-05-12 104308.jpg
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
bigben
Beiträge: 2771
Registriert: Mi Okt 12, 2016 9:09 am

Re: KPSS - Test

Beitrag von bigben »

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
---
Programmiere stets so, dass die Maxime Deines Programmierstils Grundlage allgemeiner Gesetzgebung sein könnte
Luk12Luk12
Beiträge: 25
Registriert: Di Dez 01, 2020 12:31 pm

Re: KPSS - Test

Beitrag von Luk12Luk12 »

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 :D
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
schubbiaschwilli
Beiträge: 253
Registriert: Di Jun 27, 2017 12:09 pm

Re: KPSS - Test

Beitrag von schubbiaschwilli »

Gude!

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

Dank&Gruß
Schubbiaschwilli
Luk12Luk12
Beiträge: 25
Registriert: Di Dez 01, 2020 12:31 pm

Re: KPSS - Test

Beitrag von Luk12Luk12 »

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

Danke und beste Grüße,
Lukas
bigben
Beiträge: 2771
Registriert: Mi Okt 12, 2016 9:09 am

Re: KPSS - Test

Beitrag von bigben »

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
---
Programmiere stets so, dass die Maxime Deines Programmierstils Grundlage allgemeiner Gesetzgebung sein könnte
Luk12Luk12
Beiträge: 25
Registriert: Di Dez 01, 2020 12:31 pm

Re: KPSS - Test

Beitrag von Luk12Luk12 »

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
schubbiaschwilli
Beiträge: 253
Registriert: Di Jun 27, 2017 12:09 pm

Re: KPSS - Test

Beitrag von schubbiaschwilli »

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>
> 
Luk12Luk12
Beiträge: 25
Registriert: Di Dez 01, 2020 12:31 pm

Re: KPSS - Test

Beitrag von Luk12Luk12 »

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 :lol:

Danke und beste Grüße,
Lukas
bigben
Beiträge: 2771
Registriert: Mi Okt 12, 2016 9:09 am

Re: KPSS - Test

Beitrag von bigben »

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
---
Programmiere stets so, dass die Maxime Deines Programmierstils Grundlage allgemeiner Gesetzgebung sein könnte
Antworten