[PYTHON] Dicky Fuller Test von Statistikmodellen

Dieser Artikel wurde von dem, was ich vor ein paar Tagen gepostet habe, zurückgezogen und erneut verbessert.

Es ist wichtig zu wissen, ob die zu analysierenden Zeitreihendaten einem zufälligen Spaziergang von den folgenden zwei Punkten folgen.

  1. Ist der aufgetretene Trend ein gefälschter Trend oder tritt der Trend auf?
  2. Wenn die beiden Zeitreihen eine Korrelation zeigen, handelt es sich um eine gefälschte Korrelation oder treten sie zufällig auf? Das ist der Punkt.

Viele wirtschaftliche Zeitreihen und Lagerbestände (Originalserien) folgen einem zufälligen Spaziergang. Eine Zeitreihe, die einem zufälligen Spaziergang folgt, ist eine instationäre Zeitreihe, in der sich Mittelwert und Varianz im Zeitverlauf von den vorherigen unterscheiden. Bei einem zufälligen Spaziergang ändert sich der Mittelwert wahrscheinlich über die Zeit und die Varianz nimmt proportional zur Quadratwurzel des Zeitablaufs zu. Es gibt auch verschiedene instationäre Zeitreihen. Eine davon ist die Trendstabilität, bei der der Durchschnittswert im Laufe der Zeit stetig steigt und fällt. Um festzustellen, ob es sich um einen zufälligen Spaziergang handelt, bestimmen Sie, ob die Zeitreihe eine Einheitswurzel hat. Es gibt einen Dicky-Fuller-Test als Methode. Es wird jedoch häufig darauf hingewiesen, dass die Erfassungsleistung gering ist. Also werde ich mich ein wenig darum kümmern.

Bevor Sie mit dem DF-Test fortfahren

Erstens ein Selbstrückgabemodell erster Ordnung (AR (1)) y_t=\rho y_{t-1}+u_t Wird besorgt. Wenn $ \ rho $ 1 ist, hat dieser AR (1) eine Einheitswurzel. $ y_t $ ist ein zufälliger Spaziergang. Wenn $ \ Delta y_t = y_t-y_ {t-1} $ ein stetiger Prozess mit konstantem Mittelwert und zeitlicher Varianz ist, spricht man von einer Summe erster Ordnung. Dieser Unit-Root-Prozess y_1=y_0+e_1 y_2=y_1+e_2 y_3=y_2+e_3 ​... Und wenn Sie t = 2 und t = 3 umschreiben y_2=(y_0+e_1)+e2=x_0+e1+e2 y_3=(y_0+e_1+e_2)+e_3=y_0+e_1+e_2+e_3 Es wird sein. Verallgemeinern y_t=y_0+\sum_{i=1}^te_i Ist die Summe aus Anfangswert und Zufallszahl.

Da die Zeitreihen in $ \ rho> 1 $ abweichen, setzen Sie $ \ rho <1 $. Dann macht das AR (1) -Modell dasselbe. y_1=ρ\cdot y_0+e_1 y_2=ρ\cdot y_1+e_2 y_3=ρ\cdot y_2+e_3 Umschreiben von $ y_2 $ und $ y_3 $ y_2=ρ\cdot (ρ\cdot y_0+e_1)+e2=ρ^2\cdot y_0+ρ\cdot e_1+e2 y_3=ρ\cdot (ρ^2\cdot y_0 +ρ\cdot e_1+e_2)+e_3 =ρ^3y_0+ρ^2\cdot (e_1)+ρ\cdot (e_2)+e_3 =ρ^3y_0+\sum_{i=1}^3ρ^{3-i}\cdot (e_i) Es wird sein. Verallgemeinern y_t=ρ^ty_0+\sum_{i=1}^tρ^{t-i}\cdot (e_i) Es wird sein.

Lassen Sie uns bis zu diesem Punkt mit Python überprüfen.

import numpy as np
import matplotlib.pyplot as plt
def generate(y0,rho,sigma,nsample):
    e = np.random.normal(size=nsample)
    y=[]
    y.append(rho*y0+e[0])
    for i,u in enumerate(e[1:]):
        y.append(rho*y[i]+u)
    plt.plot(y)
generate(1,1,1,100)

​image.png

Als nächstes setzen wir $ rho = 0.9 $.

generate(1,0.9,1,100)

image.png Sie können sehen, dass es zu einem bestimmten Wert konvergiert. Als nächstes setzen wir Sigma = 0 für beide.

generate(1,1,0,100)

image.png

generate(1,0.9,0,100)

image.png

Es bleibt 1 im Random Walk und konvergiert in AR (1) gegen Null. Dies ist unpraktisch. Schreiben Sie das AR (1) -Modell wie folgt um. y_t=c+\rho y_{t-1}+u_t Wird besorgt. Fügen Sie als Nächstes diese Serie hinzu. \sum_{t=2}^n y_t=\sum_{t=2}^n c+\sum_{t=2}^n\rho y_{t-1}+\sum_{t=2}^nu_t

Angenommen, der Durchschnittswert von $ y_t $ ist $ \ mu $. Wenn n groß genug ist, dann ist $ \ sum_ {t = 2} ^ nu_t = 0 $ Es wird $ \ mu = c + \ rho \ mu $. Deshalb c=\mu(1-\rho) Es wird sein. Lassen Sie uns die damit generierte Funktion neu schreiben.

def generate(y0,rho,c,sigma,nsample):
    e = np.random.normal(size=nsample)
    y=[]
    y.append(rho*y0+e[0]*sigma+c)
    for i,u in enumerate(e[1:]):
        y.append(rho*y[i]+u*sigma+c)
    plt.plot(y)
generate(1,0.9,0.1,0,100)

image.png

ist geworden. Der Durchschnittswert von AR (1) ist jetzt über die Zeit konstant. Lassen Sie uns Sigma kleiner machen und es überprüfen.

generate(1,0.9,0.1,0.01,100)

image.png

Irgendwie kann man sehen, dass es in die Mitte zurückkehrt. Das AR (1) -Modell wird stabil (schwach stabil), wenn bestimmte Bedingungen erfüllt sind. Dieses c wird manchmal als Driftrate bezeichnet. Die Driftrate ist die Rate, mit der sich der Durchschnitt ändert. Aber das ändert sich hier nicht. C ist jedoch Wenn es $ c> \ mu (1- \ rho) $ oder $ c <\ mu (1- \ rho) $ ist, ändert sich die Rate.

In der stochastischen Theorie ist die stochastische Drift die Änderung des Durchschnittswerts des stochastischen Prozesses. Es gibt eine Driftrate mit einem ähnlichen Konzept, die sich jedoch wahrscheinlich nicht ändert. Es ist ein nicht probabilistischer Parameter.

In Langzeitstudien zu Langzeitereignissen werden Zeitreihenmerkmale häufig durch Aufteilung in Trend-, periodische und probabilistische Komponenten (probabilistische Drifts) konzeptualisiert. Die Komponenten der periodischen und probabilistischen Drift werden durch Autokorrelationsanalyse und den Unterschied in den Trends identifiziert. Die Autokorrelationsanalyse hilft dabei, die richtige Phase des angepassten Modells zu identifizieren, und das Wiederholen der Differenz wandelt die stochastische Driftkomponente in weißes Rauschen um.

Trend stetig

Der trendstetive Prozess $ y_t $ folgt $ y_t = f (t) + u_t $. Wobei t die Zeit ist, f eine deterministische Funktion ist und $ u_t $ eine stationäre Wahrscheinlichkeitsvariable mit einem Langzeitmittelwert von Null ist. In diesem Fall ist der stochastische Term stationär und daher gibt es keine stochastische Drift. Die deterministische Komponente f (t) hat keinen Durchschnittswert, der langfristig konstant bleibt, aber die Zeitreihen selbst können driften. Diese deterministische Drift kann berechnet werden, indem $ y_t $ bei t zurückgegeben und aus den Daten entfernt wird, wenn ein stetiger Rest erhalten wird. Der einfachste trendsteady Prozess ist $ y_t = d \ cdot t + u_t $. $ d $ ist eine Konstante.

Gestaffelt

Der Einheitswurzelprozess (stationär) folgt $ y_t = y_ {t-1} + c + u_t $. $ u_t $ ist eine stationäre Wahrscheinlichkeitsvariable mit einem Durchschnittswert, der auf lange Sicht Null sein wird. Wobei c ein nicht-probabilistischer Driftparameter ist. Auch ohne den zufälligen Schock $ u_t $ ändert sich der Durchschnitt von y mit der Zeit um c. Diese Nichtstationarität kann aus den Daten entfernt werden, indem eine Differenz erster Ordnung genommen wird. Die abgestufte Variable $ z_t = y_t-y_ {t-1} = c + u_t $ hat einen langfristigen konstanten Mittelwert von c, so dass es keine Abweichung in diesem Mittelwert gibt. Als nächstes betrachte man c = 0. Dieser Einheitswurzelprozess ist $ y_t = y_ {t-1} + c + u_t $. Obwohl $ u_t $ ein stationärer Prozess ist, verursacht das Vorhandensein seines probabilistischen Schocks $ u_t $ eine Drift, insbesondere eine probabilistische Drift, in $ y_t $. Sobald ein Wert ungleich Null u auftritt, wird er für denselben Zeitraum in y aufgenommen. Dies ist eine Periode hinter y nach einer Periode und wirkt sich auf den y-Wert für die nächste Periode aus. y wirkt sich nacheinander auf neue y-Werte aus. Nachdem der erste Schock y beeinflusst hat, wird sein Wert permanent in den Durchschnitt von y einbezogen, was zu einer probabilistischen Drift führt. Diese Drift kann also auch beseitigt werden, indem zuerst y differenziert wird, um z zu erhalten, das nicht driftet. Diese stochastische Drift tritt auch in $ c \ ne0 $ auf. So passiert es auch mit $ y_t = y_ {t-1} + c + u_t $. Ich erklärte, dass c ein nicht-probabilistischer Driftparameter ist. Warum ist es keine eindeutige Abweichung? Tatsächlich ist die Drift, die sich aus $ c + u_t $ ergibt, Teil des stochastischen Prozesses.

numpy.random.normal Aus der numpy Referenz

numpy.random.normal(loc=0.0, scale=1.0, size=None)

loc float oder array_like von float: Mittelwert der Verteilung

scalefloat oder array_like von Floats: Standardabweichung der Verteilung

ist. Betrachten Sie nun die Bedeutung von loc. Der Durchschnittswert der aus dieser Verteilung erzeugten Proben folgt eng $ N (μ, σ ^ 2 / n) $. Mal sehen, wie es funktioniert. Wagen Sie es, loc = 1, scale = 1.0 zu setzen. Ich habe versucht, den Stichprobenmittelwert und die Standardabweichung mit n = 2 zu ermitteln. Es gibt jedoch zwei Berechnungsmethoden. Eins ist Finden Sie den Mittelwert und die Standardabweichung von $ z_t $. Der Andere ist Subtrahieren Sie 1 von jeder Beobachtung $ z_t $, um den Mittelwert und die Standardabweichung der Probe zu ermitteln. Addieren Sie dann 1 zum Stichprobenmittelwert.

import numpy as np
import matplotlib.pyplot as plt

s = np.random.normal(1, 1, 2)
print(np.mean(s),np.std(s))
o=np.ones(len(s))
Print(np.mean(s-o)+1,np.std(s-o))

# 0.9255104221128653 0.5256849930003175
# 0.9255104221128653 0.5256849930003175 #1 wurde von der Probe abgezogen und 1 zum Durchschnittswert addiert.

Das Ergebnis ist das gleiche. Es ist ersichtlich, dass np.random.normal (1,1, n) Zufallszahlen mit 1 / n + np.random.normal (0,1, n) erzeugt. Aber das ist eine Geschichte im Computer. In Wirklichkeit sollte es nicht zu unterscheiden sein, welcher die stochastische Drift erzeugt. Lassen Sie uns also 1000000 Mal Daten für n = 2 und n = 100 generieren und sehen, wie das Ergebnis aussieht.

Verteilung der arithmetischen Mittelwerte von Stichproben gestaffelter stetiger Prozesse

Die abgestufte Variable $ z_t = y_t-y_ {t-1} = c + u_t $ hat einen langfristigen konstanten Mittelwert von c. Der Stichprobenmittelwert folgt genau $ N (μ, σ ^ 2 / n) $, wenn n groß ist. Daher ist der Durchschnittswert nicht konstant. In diesem Fall gibt es jedoch keine Abweichung im Durchschnittswert. Es hat einen Durchschnittswert, der langfristig konstant ist.

m1=[]
trial=1000000
for i in range(trial):
    s = np.random.normal(1, 1, 2)
    m1.append(np.mean(s))
m2=[]
for i in range(trial):
    s = np.random.normal(1, 1, 100)
    m2.append(np.mean(s))
bins=np.linspace(np.min(s),np.max(s),1000)
plt.hist(m1,bins=bins,label='n=2')
plt.hist(m2,bins=bins,label='n=100',alpha=0.5)
plt.legend()
plt.show()

image.png

Sowohl der Stichprobenmittelwert als auch die Standardabweichung bilden eine Verteilung und werden nicht eindeutig bestimmt. N muss viel größer sein, um eindeutig zu sein. Daher werden die Merkmale der Zeitreihen in Trendkomponenten, periodische Komponenten und probabilistische Komponenten (probabilistische Drift) zerlegt und konzeptualisiert. Beim Umgang mit probabilistischer Drift und Driftrate ist jedoch Vorsicht geboten.

Probabilistische Drift tritt im Einheitswurzelprozess auf. Selbst wenn es eine bestimmte Driftrate gibt, ist sie nicht zu unterscheiden.

Verteilung von y mit Einheitswurzeln am letzten Tag

Im Einheitswurzelprozess nimmt die Standardabweichung mit zunehmender Anzahl von Proben zu. Der Durchschnittswert von y mit einer Stichprobengröße von n ist Null bei c = 0, aber die Verteilung erweitert sich, wenn die Standardabweichung $ \ sqrt {n \ cdot \ sigma ^ 2} $ zunimmt.

m1=[]
trial=1000000
for i in range(trial):
    s = np.cumsum(np.random.normal(0, 1, 2))
    m1.append(s[-1])
m2=[]
for i in range(trial):
    s = np.cumsum(np.random.normal(0, 1, 100))
    m2.append(s[-1])
bins=np.linspace(np.min(m2),np.max(m2),1000)
plt.hist(m1,bins=bins,label='n=2')
plt.hist(m2,bins=bins,label='n=100',alpha=0.5)
plt.legend()
plt.show()

image.png

Die Verteilung unterscheidet sich von der vorherigen. Je größer die Stichprobe ist, desto breiter ist die Verteilung. Versuchen wir, die eigene Serie zu zeichnen.

plt.plot(s)

image.png

In diesem Fall ändert sich der Pegel jeder Probe schrittweise. Diese Änderung ist auf einen stochastischen Schock $ u_t $ zurückzuführen. Was hinzugefügt wurde, wird aufgezeichnet. Eine Reihe von stochastischen Abweichungen ist dargestellt. Dies kann einen großen Trend bilden. Sie entspricht der Fläche an beiden Enden in der Verteilungskarte von n = 100. Welcher Bereich von der Definition des Trends abhängt. Dieser Trend ist ein probabilistischer Trend. Dies unterscheidet sich deutlich von den deterministischen Trends, die durch Trendstabilität entstehen.

Selbstrückgabeprozess erster Ordnung

Der Selbstrückgabeprozess erster Ordnung kann wie folgt umgeschrieben werden. x_1=ρ\cdot x_0+y_1+e_1 x_2=ρ\cdot x_1+y_2+e_2 x_3=ρ\cdot x_2+y_3+e_3

x_2=ρ\cdot (ρ\cdot x_0+y_1+e_1)+e2=ρ^2\cdot x_0+ρ\cdot (y_1+e1)+y_2+e2 x_3=ρ\cdot (ρ\cdot ρ\cdot x_0+ρ\cdot (y_1+e_1)+y_2+e_2)+y_3+e_3 =ρ^3x_0+ρ^2\cdot (y_1+e_1)+ρ\cdot (y_2+e_2)+y_3+e_3 =ρ^3x_0+\sum_{i=1}^n\rho^{n-i}(y_i+e_i) Wobei $ y_i = c + d \ cdot i $

Wenn $ \ rho = 1 $ ist, handelt es sich um einen zufälligen Walk-Prozess.

#Initialisieren
import numpy as np
import statsmodels.api as sm
import matplotlib.pyplot as plt
from scipy.stats import norm
from statsmodels.tsa.stattools import adfuller
#np.random.seed(0)

# AR(1)Prozessgenerierung
#x0 ist der Anfangspreis
#nsample ist die Anzahl der Perioden in einem Versuch
# rho=1: random walk
## c=c1-Rho ist Drift, c1=Wenn es 1 ist, gibt es keine Drift. Random Walk d ist ein eindeutiger Trend
# rho<Selbstrückgabemodell 1 bis 1. Ordnung
##Die Auswirkungen von x0, c und d verschwinden mit der Zeit
def ar1(rho,x0,c1,d,sigma,nsample):
    e = np.random.normal(size=nsample)#random variable
    rhoo=np.array([rho ** num for num in range(nsample)]) #effects of rho
    rhoo_0=np.array([rho**(num+1) for num in range(nsample)]) #effects of rho
    t=np.arange(nsample) #deterministic trend
    one=np.ones(nsample) #effect of x0
    x00=x0*one*rhoo_0 #effect of rho to x0
    c=c1-rho # constant of autoregressive mdoel
    c0=np.cumsum(c*one*rhoo) #effect of rho to c
    d0=(d*t*rhoo) # effect of rho to determinist trend
    y=np.cumsum(rhoo*e*sigma)+x00+c0+d0
    return y

Verteilung von y am letzten Tag des AR (1) -Prozesses

m1=[]
trial=1000000
nsample=250
x0=1
for i in range(trial):
    s = np.cumsum(np.random.normal(0, 1, nsample))
    m1.append(s[-1]+x0)
def generate(rho,x0,c1,sigma,nsample,trial):
    m=[]
    for i in range(trial):
        e = np.random.normal(0, 1, nsample)
        rhoo=np.array([rho ** num for num in range(nsample)]) #effects of rho
        rhoo_0=np.array([rho**(num+1) for num in range(nsample)]) #effects of rho
        one=np.ones(nsample) #effect of x0
        x00=x0*one*rhoo_0 #effect of rho to x0
        c=c1-rho # constant of autoregressive mdoel
        c0=np.cumsum(c*one*rhoo) #effect of rho to c
        s=np.cumsum(rhoo*e*sigma)+x00+c0
        m.append(s[-1])
    return m
c1=1
x0=1
rho=0.99
sigma=1
m2=generate(rho,x0,c1,sigma,nsample,trial)
rho=0.9
m3=generate(rho,x0,c1,sigma,nsample,trial)
rho=0.5
m4=generate(rho,x0,c1,sigma,nsample,trial)
bins=np.linspace(np.min(m1),np.max(m1),1000)
plt.hist(m1,bins=bins,label='rw')
plt.hist(m2,bins=bins,label='ar1 rho=0.99',alpha=0.5)
plt.hist(m3,bins=bins,label='ar1 rho=0.9',alpha=0.5)
plt.hist(m4,bins=bins,label='ar1 rho=0.5',alpha=0.5)
plt.legend()
plt.show()

image.png Im Self-Return-Modell erster Ordnung können Sie sehen, dass $ \ rho $, wenn es kleiner als 1 wird, um den Anfangswert herum wandert. Wenn sich $ \ rho $ 1 nähert, wird es schwierig, es von einem zufälligen Spaziergang zu unterscheiden. Die folgende Abbildung gilt für n = 7. In ungefähr 7 Tagen wird die Streuung von $ \ rho $ ähnlich sein.

image.png

Ermitteln der Leistung (Einführung in das Testen statistischer Hypothesen anhand von Statistikmodellen)

Im Hypothesentest werden die Nullhypothese und die Alternativhypothese erstellt. Die Nullhypothese besteht aus einer Population. Da die Bevölkerung im Allgemeinen unbekannt ist, kann eine Schätzung verwendet werden. Die Proben werden als aus dieser Population extrahiert behandelt. In der folgenden Abbildung ist die Verteilung der Bevölkerung durch die blau gepunktete Linie dargestellt. Die Stichprobe gibt die Stichprobenverteilung dieser Statistik an. Es wird durch eine orange Linie angezeigt. Vergleichen Sie die beiden. Wenn die beiden Verteilungen nahe beieinander liegen, ist es schwierig zu beurteilen, dass die beiden Verteilungen unterschiedlich sind, und wenn sie weit voneinander entfernt sind, kann beurteilt werden, dass sie unterschiedlich sind. image.png Diese Beurteilung kann anhand des Signifikanzniveaus erfolgen. Das Signifikanzniveau $ \ alpha $ ist der Bereich ganz rechts in der Bevölkerung. 10%, 5% usw. werden verwendet. Der Bereich rechts von der blauen vertikalen Linie in der blauen Verteilung. Die blaue vertikale Linie zeigt die Grenze, an der die Nullhypothese verworfen wird. Der Bereich innerhalb dieser Linie wird als Akzeptanzbereich bezeichnet, und der Bereich außerhalb dieser Linie wird als Ablehnungsbereich bezeichnet.

Sobald die Stichprobe erhalten wurde, können Statistiken wie Mittelwert und Varianz berechnet werden. Der p-Wert ist die Wahrscheinlichkeit der Zurückweisung, wenn die Nullhypothese mit der Statistik als Untergrenze zurückgewiesen wird, die aus der Verteilung der Population berechnet wird. Als nächstes ist es in der Abbildung der Bereich rechts von der roten vertikalen Linie. Wenn dieser Wert kleiner als das Signifikanzniveau ist, das das Kriterium ist, wird beurteilt, dass die Stichprobe (Statistik) aus einer Verteilung erhalten wurde, die sich von der Bevölkerung (Bevölkerung) unterscheidet. In der Abbildung befindet sich die rote vertikale Linie links von der blauen vertikalen Linie. Wenn jedoch der Durchschnittswert der orangefarbenen Stichprobenverteilung rechts von der blauen vertikalen Linie liegt, ist klar, dass ihr p-Wert unter dem Signifikanzniveau liegt. image.png Die Fähigkeit, eine falsche Nullhypothese korrekt abzulehnen, wird als Erkennung bezeichnet. Das Erkennen der Leistung $ \ beta $ entspricht dem Bereich rechts von der orangefarbenen Verteilung der blauen vertikalen Linien. Wenn der $ \ beta $ -Teil groß ist, steigt die Wahrscheinlichkeit, dass die alternative Hypothese korrekt ist, wenn die Nullhypothese verworfen wird. image.png Um dieses $ \ beta $ zu erhöhen, ist es möglich, dass die beiden Verteilungen weit voneinander entfernt sind oder die Positionen der Verteilungszentren erheblich voneinander abweichen. Das andere ist die Größe der Probe. Dann wird die Form der Verteilung verengt, und wenn die Größe der Probe groß genug ist, um die Verteilung als Linie zu betrachten, werden die beiden geraden Linien verglichen. image.png

Dicky-Fuller-Test

Ein einfacher AR (1) -Prozess $ y_t= \rho y_{t-1}+u_t$ Und subtrahiere $ y_ {t-1} $ von beiden Seiten y_t-y_{t-1}=(\rho-1)y_{t-1}+u_t Wird besorgt. Als $ \ delta = \ rho-1 $ lautet die Nullhypothese, dass das AR-Modell eine Einheitswurzel enthält. H_0:\delta=0 H_1: \delta<0 Es wird sein.

Es gibt drei Arten.

  1. Zufälliges Gehen ohne Drift \Delta y_t=\delta y_{t-1}+u_t
  2. Zufälliger Spaziergang mit Drift \Delta y_t=c+\delta y_{t-1}+u_t
  3. Zufälliger Spaziergang mit Drift + Zeittrend \Delta y_t=c+c\cdot t+\delta y_{t-1}+u_t Dies sind Unterschiede erster Ordnung auf der linken Seite, was bedeutet, dass die Unterschiede Abweichungen oder Zeittrends aufweisen.
# Dicky-Voller Test
def adfcheck(a0,mu,sigma,nsample,trial):
    prw=0
    prwc=0
    prwct=0
    for i in range(trial):
        y=ar1(a0,1,mu,sigma,nsample)
        rw=sm.tsa.adfuller(y,maxlag=0,regression='nc')[1]
        rwc=sm.tsa.adfuller(y,maxlag=0,regression='c')[1]
        rwct=sm.tsa.adfuller(y,maxlag=0,regression='ct')[1]
        if rw>0.05:
            prw+=1
        if rwc>0.05:
            prwc+=1
        if rwct>0.05:
            prwct+=1
    print('nsample',nsample,prw/trial,prwc/trial,prwct/trial)

Landau Walk Test

Zufälliger Spaziergang ohne Drift

rho=1
c1=1
c1=(c1-1)/250+1
d=0
d=d/250
sigma=0.4/np.sqrt(250)
trial=10000    
for i in [7,20,250]:
    adfcheck(rho,x0,c1,d,sigma,i,trial)

# nsample 7 0.9448 0.8878 0.8528
# nsample 20 0.9635 0.9309 0.921
# nsample 250 0.9562 0.9509 0.9485

Zufälliger Spaziergang mit Drift

rho=1
c1=1.1
c1=(c1-1)/250+1
d=0
d=d/250
sigma=0.4/np.sqrt(250)
trial=10000    
for i in [7,20,250]:
    adfcheck(rho,x0,c1,d,sigma,i,trial)

# nsample 7 0.9489 0.8889 0.8524
# nsample 20 0.9649 0.9312 0.921
# nsample 250 0.9727 0.9447 0.9454

Random Walk + Drift + deterministischer Trend

rho=1
c1=1.1
c1=(c1-1)/250+1
d=0.9
d=d/250
sigma=0.4/np.sqrt(250)
trial=10000    
for i in [7,20,250]:
    adfcheck(rho,x0,c1,d,sigma,i,trial)

# nsample 7 0.9971 0.948 0.8581
# nsample 20 1.0 1.0 0.983
# nsample 250 1.0 1.0 0.9998

Wenn der zufällige Spaziergang abgelehnt werden soll

AR (1) Modell ohne Drift

rho=0.99
c1=rho
#c1=(c1-1)/250+1
d=0
d=d/250
sigma=0.4/np.sqrt(250)
trial=10000    
for i in [7,20,250]:
    adfcheck(rho,x0,c1,d,sigma,i,trial)

# nsample 7 0.7797 0.9012 0.8573
# nsample 20 0.6188 0.9473 0.9225
# nsample 250 0.0144 0.7876 0.944

AR (1) Modell mit Drift

rho=0.90
c1=1.1
c1=(c1-1)/250+1
d=0
d=d/250
sigma=0.4/np.sqrt(250)
trial=10000    
for i in [7,20,250]:
    adfcheck(rho,x0,c1,d,sigma,i,trial)

# nsample 7 0.9715 0.8868 0.8539
# nsample 20 0.9989 0.9123 0.912
# nsample 250 1.0 0.0306 0.1622

AR (1) Modell + Drift + deterministischer Trend

rho=0.90
c1=1.1
c1=(c1-1)/250+1
d=1
d=d/250
sigma=0.4/np.sqrt(250)
trial=10000    
for i in [7,20,250]:
    adfcheck(rho,x0,c1,d,sigma,i,trial)

# nsample 7 0.999 0.9404 0.8576
# nsample 20 1.0 1.0 0.918
# nsample 250 1.0 1.0 0.0044

Wir können sehen, dass wir ein Modell mit den richtigen Spezifikationen verwenden sollten, dass die Stichprobengröße groß sein sollte und dass ρ weit von 1 entfernt sein sollte, wenn die Nullhypothese verworfen wird.

Über den Einfluss von Lärm

Verwenden wir die Cauchy-Verteilung, um Fettschwanzgeräusche zu erzeugen und um zu sehen, wie sie sich auf die Zeitreihen auswirken. Schauen wir uns zunächst die Verteilung der Cauchy-Verteilung an.

from scipy.stats import cauchy
trial=1000000
m1=[]
for i in range(trial):
    s = np.random.normal(0, 1, 250)
    m1.append(np.mean(s))
m2=[]
for i in range(trial):
    s=[]
    i=1
    while i<=250:
        s0 = (cauchy.rvs(0,1,1))
        if abs(s0)<10:
            s.append(s0)
            i+=1
    m2.append(np.mean(s))
bins=np.linspace(np.min(m2),np.max(m2),1000)
plt.hist(m1,bins=bins,label='normal')
plt.hist(m2,bins=bins,label='abs(cauchy) <10',alpha=0.5)
plt.legend()
plt.show()

image.png

Sie können sehen, dass der fette Schwanz stark ist. Betrachten Sie dies als Lärm.

Bilden Sie als Nächstes die Summe der Zufallszahlen aus der Koshi-Verteilung und sehen Sie sich die Verteilung des letzten Werts an.

from scipy.stats import cauchy
m1=[]
trial=1000000
for i in range(trial):
    s = np.cumsum(np.random.normal(0, 1, 250))
    m1.append(s[-1])
m2=[]
for i in range(trial):
    s=[]
    i=1
    while i<=250:
        s0 = (cauchy.rvs(0,1,1))
        if abs(s0)<10:
            s.append(s0)
            i+=1
    m2.append(np.cumsum(s)[-1])
bins=np.linspace(np.min(m2),np.max(m2),1000)
plt.hist(m1,bins=bins,label='normal')
plt.hist(m2,bins=bins,label='cauchy',alpha=0.5)
plt.legend()
plt.show()

image.png Immerhin können Sie sehen, dass der Spread größer als die Normalverteilung ist.

Als nächstes erstellen wir eine AR1-Funktion, die die Cauchy-Verteilung verwendet.

def ar1_c(rho,x0,c1,d,sigma,nsample):
    e = cauchy.rvs(loc=0,scale=1,size=nsample)    
    c=c1-rho # constant of autoregressive mdoel
    y=[]
    y.append(e[0]+x0*rho+c+d)
    for i,ee in enumerate(e[1:]):
        y.append(rho*y[i]+ee+c+d*(i+1))
    return y

Der DF-Test wird unter Verwendung der erstellten Selbstrückgabezeitreihen erster Ordnung durchgeführt.

def adfcheck_c(rho,x0,c1,d,sigma,nsample,trial):
    prw=0
    prwc=0
    prwct=0
    for i in range(trial):
        y=ar1_c(rho,x0,c1,d,sigma,nsample)
        rw=sm.tsa.adfuller(y,maxlag=0,regression='nc')[1]
        rwc=sm.tsa.adfuller(y,maxlag=0,regression='c')[1]
        rwct=sm.tsa.adfuller(y,maxlag=0,regression='ct')[1]
        if rw>0.05:
            prw+=1
        if rwc>0.05:
            prwc+=1
        if rwct>0.05:
            prwct+=1
    print('nsample',nsample,prw/trial,prwc/trial,prwct/trial)

rho=0.90
c1=1.1
c1=(c1-1)/250+1
d=1
d=d/250
sigma=0.4/np.sqrt(250)
trial=10000    
for i in [7,20,250]:
    adfcheck_c(rho,x0,c1,d,sigma,i,trial)

# nsample 7 0.973 0.8775 0.8377
# nsample 20 0.9873 0.9588 0.9143
# nsample 250 0.8816 0.8425 0.0897

Aufgrund des Einflusses des Fettschwanzes handelt es sich in diesem Fall um einen Zeittrend mit Drift. Wenn die Stichprobengröße jedoch 250 beträgt, kann die Nullhypothese bei einem Signifikanzniveau von 5% nicht verworfen werden. Wenn rho = 0,85 ist, ist das Ergebnis wie folgt.

nsample 7 0.9717 0.8801 0.8499 nsample 20 0.987 0.9507 0.903 nsample 250 0.8203 0.7057 0.0079

Wenn die Stichprobengröße 250 beträgt, kann die Nullhypothese bei einem Signifikanzniveau von 5% verworfen werden.

Als nächstes schneiden wir den fetten Schwanzteil der Cauchy-Verteilung ab.

def cau(cut,nsample):
    s=[]
    i=1
    while i<=nsample:
        s0 = cauchy.rvs(0,1,1)
        if abs(s0)<cut:
            s.append(s0[0])
            i+=1
    return s

def ar1_c2(rho,x0,c1,d,sigma,nsample):
    e = cau(10,nsample)
    c=c1-rho # constant of autoregressive mdoel  
    y=[]
    y.append(e[0]*sigma+x0*rho+c+d)
    for i,ee in enumerate(e[1:]):
        y.append(rho*y[i]+ee*sigma+c+d*(i+1))
    return y

def adfcheck_c2(rho,x0,c1,d,sigma,nsample,trial):
    prw=0
    prwc=0
    prwct=0
    for i in range(trial):
        y=ar1_c2(rho,x0,c1,d,sigma,nsample)
        rw=sm.tsa.adfuller(y,maxlag=0,regression='nc')[1]
        rwc=sm.tsa.adfuller(y,maxlag=0,regression='c')[1]
        rwct=sm.tsa.adfuller(y,maxlag=0,regression='ct')[1]
        if rw>0.05:
            prw+=1
        if rwc>0.05:
            prwc+=1
        if rwct>0.05:
            prwct+=1
    print('nsample',nsample,prw/trial,prwc/trial,prwct/trial)

rho=0.90
c1=1.1
c1=(c1-1)/250+1
d=1
d=d/250
sigma=0.4/np.sqrt(250)
trial=10000    
for i in [7,20,250]:
    adfcheck_c2(rho,x0,c1,d,sigma,i,trial)

# nsample 7 0.9969 0.9005 0.8462
# nsample 20 1.0 0.9912 0.9182
# nsample 250 1.0 1.0 0.0931

Das Ergebnis war anders als das vorherige. Nur Zeittrends mit Drift können die Nullhypothese ablehnen, aber das Signifikanzniveau beträgt 10%. Wenn rho = 0,85

nsample 7 0.9967 0.9012 0.8466 nsample 20 1.0 0.9833 0.9121 nsample 250 1.0 1.0 0.0028

Das Ergebnis ist.

Versuchen Sie als nächstes den erweiterten Dicky Fuller-Test.

\Delta y_t=c+c\cdot t+\delta y_{t-1}+\delta_1\Delta y_{t-1}+,\cdots,+\delta_{p-1}\Delta y_{t-p+1}+u_t Durch Erhöhen der Differenz erster Ordnung der Verzögerung werden strukturelle Effekte wie Autokorrelation beseitigt.

def adfcheck_c3(rho,x0,c1,d,sigma,nsample,trial):
    prw=0
    prwc=0
    prwct=0
    for i in range(trial):
        y=ar1_c2(rho,x0,c1,d,sigma,nsample)
        rw=sm.tsa.adfuller(y,regression='nc')[1]
        rwc=sm.tsa.adfuller(y,regression='c')[1]
        rwct=sm.tsa.adfuller(y,regression='ct')[1]
        if rw>0.05:
            prw+=1
        if rwc>0.05:
            prwc+=1
        if rwct>0.05:
            prwct+=1
    print('nsample',nsample,prw/trial,prwc/trial,prwct/trial)

rho=0.90
c1=1.1
c1=(c1-1)/250+1
d=1
d=d/250
sigma=0.4/np.sqrt(250)
trial=10000    
for i in [7,20,250]:
    adfcheck_c3(rho,x0,c1,d,sigma,i,trial)

# nsample 7 0.94 0.8187 0.8469
# nsample 20 1.0 0.873 0.7969
# nsample 250 1.0 1.0 0.1253

Es kann nicht richtig mit rho = 0,90 beurteilt werden. Wenn rho = 0,85 ist, ist das Ergebnis

nsample 7 0.9321 0.8109 0.8429 nsample 20 1.0 0.8512 0.7942 nsample 250 1.0 1.0 0.0252

Ich konnte den Drift + Zeit-Trend ablehnen.

Erkennung von Drift und Zeittrend, wenn die Standardabweichung Null oder klein ist

#Nur Drift, Standardabweichung ist Null
rho=0
c1=1.1
c1=(c1-1)/250+1
d=0
d=d/250
sigma=0#.4/np.sqrt(250)
trial=1#0000    
for i in [7,20,250]:
    adfcheck(rho,x0,c1,d,sigma,i,trial)

nsample 7 0.0 0.0 0.0 nsample 20 0.0 0.0 0.0 nsample 250 0.0 0.0 0.0

Die Nullhypothese wird überhaupt abgelehnt, wenn es nur eine Drift gibt.

#Drift und Zeittrend, Standardabweichung ist Null
rho=0
c1=1.1
c1=(c1-1)/250+1
d=1
d=d/250
sigma=0#.4/np.sqrt(250)
trial=1#0000    
for i in [7,20,250]:
    adfcheck(rho,x0,c1,d,sigma,i,trial)

nsample 7 1.0 1.0 0.0 nsample 20 1.0 1.0 0.0 nsample 250 1.0 1.0 0.0 AR (1) + Drift + Zeittrend lehnt die Nullhypothese unabhängig von der Stichprobengröße ab

Fügen Sie als Nächstes der Drift eine Standardabweichung hinzu.

# drift+ very small sigma
rho=0
c1=1.1
c1=(c1-1)/250+1
d=0
d=d/250
sigma=0.004#.4/np.sqrt(250)
trial=10000    
for i in [7,20,250]:
    adfcheck(rho,x0,c1,d,sigma,i,trial)

nsample 7 0.9986 0.6379 0.7223 nsample 20 1.0 0.0396 0.1234 nsample 250 1.0 0.0 0.0

Fügen Sie als Nächstes dem Trend Drift + Zeit eine Standardabweichung hinzu.

# drift+ very small sigma
rho=0
c1=1.1
c1=(c1-1)/250+1
d=1
d=d/250
sigma=0.004#.4/np.sqrt(250)
trial=10000    
for i in [7,20,250]:
    adfcheck(rho,x0,c1,d,sigma,i,trial)

nsample 7 1.0 0.9947 0.6877 nsample 20 1.0 1.0 0.1095 nsample 250 1.0 1.0 0.0

Referenz

Campbell, J. Y.; Perron, P. (1991). "Pitfalls and Opportunities: What Macroeconomists Should Know about Unit Roots"

Stock J. Unit Roots, Structural Breaks, and Trends. In: Engle R, McFadden D Handbook of Econometrics. Amsterdam: Elsevier ; 1994. pp. 2740-2843.

MacKinnon, J.G. 2010. “Critical Values for Cointegration Tests.” Queen’s University, Dept of Economics, Working Papers.

statsmodels.tsa.stattools.adfuller

Unit root

[Dicky-Fuller-Test](https://ja.wikipedia.org/wiki/%E3%83%87%E3%82%A3%E3%83%83%E3%82%AD%E3%83%BC% E2% 80% 93% E3% 83% 95% E3% 83% A9% E3% 83% BC% E6% A4% 9C% E5% AE% 9A)

Stochastic drift

http://web.econ.ku.dk/metrics/Econometrics2_05_II/Slides/08_unitroottests_2pp.pdf

Recommended Posts

Dicky Fuller Test von Statistikmodellen
Primzahlbeurteilung durch Python