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.
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.
Erstens ein Selbstrückgabemodell erster Ordnung (AR (1))
Da die Zeitreihen in $ \ rho> 1 $ abweichen, setzen Sie $ \ rho <1 $. Dann macht das AR (1) -Modell dasselbe.
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)
Als nächstes setzen wir $ rho = 0.9 $.
generate(1,0.9,1,100)
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)
generate(1,0.9,0,100)
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.
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
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)
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)
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.
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.
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.
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()
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.
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()
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)
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.
Der Selbstrückgabeprozess erster Ordnung kann wie folgt umgeschrieben werden.
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
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()
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.
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.
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.
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.
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.
Ein einfacher AR (1) -Prozess
$ y_t= \rho y_{t-1}+u_t$
Und subtrahiere $ y_ {t-1} $ von beiden Seiten
Es gibt drei Arten.
# 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)
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
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
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
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
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
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.
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()
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()
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.
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.
#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
statsmodels.tsa.stattools.adfuller
[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)
http://web.econ.ku.dk/metrics/Econometrics2_05_II/Slides/08_unitroottests_2pp.pdf