Test de Dicky Fuller par des modèles statistiques

Cet article a été retiré de ce que j'avais publié il y a quelques jours et a été amélioré à nouveau.

Il est important de savoir si les données de séries chronologiques à analyser suivent une marche aléatoire à partir des deux points suivants.

  1. La tendance qui s'est produite est-elle une fausse tendance ou la tendance se produit-elle?
  2. Lorsque les deux séries chronologiques montrent une corrélation, s'agit-il d'une fausse corrélation ou se produisent-elles? C'est le but.

De nombreuses séries chronologiques économiques et niveau de stock (série originale) suivent une marche aléatoire. Une série chronologique qui suit une marche aléatoire est une série chronologique non stationnaire dans laquelle la moyenne et la variance diffèrent des précédentes au fil du temps. Dans une marche aléatoire, la moyenne change de manière probabiliste avec le temps et la variance augmente proportionnellement à la racine carrée du passage du temps. Il existe également diverses séries chronologiques non stationnaires. L'un d'eux est la stabilité de la tendance dans laquelle la valeur moyenne augmente et diminue régulièrement avec le temps. Pour déterminer s'il s'agit d'une marche aléatoire, déterminez si la série chronologique a une racine unitaire. Il existe un test Dicky Fuller comme méthode. Cependant, il est souvent souligné que la puissance de détection est faible. Je vais donc l'examiner un peu.

Avant de passer au test DF

Tout d'abord, un modèle d'auto-retour du premier ordre (AR (1)) y_t=\rho y_{t-1}+u_t ça ira. Si $ \ rho $ vaut 1, cet AR (1) est dit avoir une racine unitaire. $ y_t $ est une marche aléatoire. Lorsque $ \ Delta y_t = y_t-y_ {t-1} $ est un processus constant avec une moyenne et une variance constantes dans le temps, on parle de somme du premier ordre. Ce processus de racine unitaire y_1=y_0+e_1 y_2=y_1+e_2 y_3=y_2+e_3 ​... Et si vous réécrivez t = 2 et t = 3 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 Ce sera. Généraliser y_t=y_0+\sum_{i=1}^te_i Est la somme de la valeur initiale et du nombre aléatoire.

Puisque la série temporelle diverge dans $ \ rho> 1 $, définissez $ \ rho <1 $. Ensuite, le modèle AR (1) fera de même. y_1=ρ\cdot y_0+e_1 y_2=ρ\cdot y_1+e_2 y_3=ρ\cdot y_2+e_3 Réécriture de $ y_2 $ et $ 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) Ce sera. Généraliser y_t=ρ^ty_0+\sum_{i=1}^tρ^{t-i}\cdot (e_i) Ce sera.

Vérifions jusqu'à ce point en utilisant python.

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

Ensuite, définissons $ rho = 0,9 $.

generate(1,0.9,1,100)    

image.png Vous pouvez voir qu'il converge vers une certaine valeur. Ensuite, définissons sigma = 0 pour les deux.

generate(1,1,0,100)    

image.png

generate(1,0.9,0,100)    

image.png

Il reste 1 sur la marche aléatoire et converge vers zéro sur AR (1). Ceci n'est pas pratique, alors réécrivez le modèle AR (1) comme suit. y_t=c+\rho y_{t-1}+u_t ça ira. Ensuite, ajoutez cette série. \sum_{t=2}^n y_t=\sum_{t=2}^n c+\sum_{t=2}^n\rho y_{t-1}+\sum_{t=2}^nu_t

En supposant que la valeur moyenne de $ y_t $ est $ \ mu $, si n est assez grand, alors $ \ sum_ {t = 2} ^ nu_t = 0 $ Cela devient $ \ mu = c + \ rho \ mu $. Donc c=\mu(1-\rho) Ce sera. Réécrivons la fonction generate en utilisant ceci.

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

est devenu. La valeur moyenne de AR (1) est désormais constante dans le temps. Rendons sigma plus petit et vérifions-le.

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

image.png

D'une manière ou d'une autre, vous pouvez voir qu'il revient au centre. Le modèle AR (1) devient stable (faible stable) si certaines conditions sont remplies. Ce c est parfois appelé le taux de dérive. La vitesse de dérive est la vitesse à laquelle la moyenne change. Mais ça ne change pas ici. Cependant, c est Si c'est $ c> \ mu (1- \ rho) $ ou $ c <\ mu (1- \ rho) $, le taux changera.

Dans la théorie stochastique, la dérive stochastique est le changement de la valeur moyenne du processus stochastique. Il existe un taux de dérive avec un concept similaire, mais cela ne change pas de manière probabiliste. C'est un paramètre non probabiliste.

Dans les études à long terme sur les événements à long terme, les caractéristiques des séries chronologiques sont souvent conceptualisées en les décomposant en composantes tendance, périodiques et probabilistes (dérives probabilistes). Les composantes de la dérive périodique et probabiliste sont identifiées par l'analyse d'autocorrélation et la différence des tendances. L'analyse d'autocorrélation permet d'identifier la phase correcte du modèle ajusté et la répétition de la différence transforme la composante de dérive stochastique en bruit blanc.

Tendance stable

Le processus de tendance stable $ y_t $ suit $ y_t = f (t) + u_t $. Où t est le temps, f est une fonction déterministe et $ u_t $ est une variable de probabilité stationnaire avec une moyenne à long terme de zéro. Dans ce cas, le terme stochastique est stationnaire et il n'y a donc pas de dérive stochastique. La composante déterministe f (t) n'a pas de valeur moyenne qui reste constante sur le long terme, mais la série temporelle elle-même peut dériver. Cette dérive déterministe peut être calculée en renvoyant $ y_t $ à t et en la supprimant des données une fois qu'un résidu stable est obtenu. Le processus de stabilisation de tendance le plus simple est $ y_t = d \ cdot t + u_t $. $ d $ est une constante.

Échelonné

Le processus de racine unitaire (stationnaire) suit $ y_t = y_ {t-1} + c + u_t $. $ u_t $ est une variable de probabilité stationnaire avec une valeur moyenne qui sera nulle à long terme. Où c est un paramètre de dérive non probabiliste. Même sans le choc aléatoire $ u_t $, la moyenne de y change de c dans le temps. Cette non-stationnarité peut être supprimée des données en prenant une différence de premier ordre. La variable graduée $ z_t = y_t-y_ {t-1} = c + u_t $ a une moyenne constante à long terme de c, il n'y a donc pas de dérive dans cette moyenne. Ensuite, considérons c = 0. Ce processus racine unitaire est $ y_t = y_ {t-1} + c + u_t $. Bien que $ u_t $ soit un processus stationnaire, la présence de son choc probabiliste $ u_t $ provoque une dérive, notamment probabiliste, dans $ y_t $. Une fois qu'une valeur u non nulle se produit, elle est incluse dans y pour la même période. Ce sera une période derrière y après une période et affectera la valeur y pour la période suivante. y affecte les nouvelles valeurs y les unes après les autres. Par conséquent, après que le premier choc affecte y, sa valeur est incorporée en permanence dans la moyenne de y, résultant en une dérive probabiliste. Donc, cette dérive peut également être éliminée en différenciant d'abord y pour obtenir z qui ne dérive pas. Cette dérive stochastique se produit également dans $ c \ ne0 $. Cela arrive donc aussi avec $ y_t = y_ {t-1} + c + u_t $. J'ai expliqué que c est un paramètre de dérive non probabiliste. Pourquoi n'est-ce pas une dérive définitive? En fait, la dérive qui résulte de $ c + u_t $ fait partie du processus stochastique.

numpy.random.normal De la référence numpy

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

loc float ou array_like of floats: valeur moyenne de la distribution

scalefloat ou array_like of floats: écart type de la distribution

est. Considérons maintenant la signification de loc. La valeur moyenne des échantillons générés à partir de cette distribution suit de près $ N (μ, σ ^ 2 / n) $. Voyons voir comment ça fonctionne. Osez définir loc = 1, échelle = 1.0. J'ai essayé de trouver la moyenne et l'écart type de l'échantillon avec n = 2. Cependant, il existe deux méthodes de calcul. L'un est Trouvez la moyenne et l'écart type de $ z_t $. L'autre est Soustrayez 1 de chaque observation $ z_t $ pour trouver la moyenne et l'écart type de l'échantillon. Ajoutez ensuite 1 à la moyenne de l'échantillon.

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 a été soustrait de l'échantillon et 1 a été ajouté à la valeur moyenne.

Le résultat est le même. On peut voir que np.random.normal (1,1, n) génère des nombres aléatoires avec 1 / n + np.random.normal (0,1, n). Mais c'est une histoire dans l'ordinateur. En réalité, il devrait être impossible de distinguer celui qui crée la dérive stochastique. Générons donc des données 1000000 fois pour n = 2 et n = 100 et voyons à quoi ressemble le résultat.

Distribution des moyennes arithmétiques d'échantillons de processus stables échelonnés

La variable graduée $ z_t = y_t-y_ {t-1} = c + u_t $ a une moyenne constante à long terme de c. La moyenne de l'échantillon suit de près $ N (μ, σ ^ 2 / n) $ si n est grand. Par conséquent, la valeur moyenne n'est pas constante. Cependant, dans ce cas, il n'y a pas de dérive de la valeur moyenne. Il a une valeur moyenne constante sur le long terme.

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

La moyenne de l'échantillon et l'écart type forment une distribution et ne sont pas déterminés de manière unique. N doit être beaucoup plus grand pour être unique. Par conséquent, les caractéristiques de la série chronologique sont décomposées en composantes de tendance, composantes périodiques et composantes probabilistes (dérive probabiliste) pour les conceptualiser, mais il faut faire attention lors de la gestion de la dérive probabiliste et du taux de dérive.

Une dérive probabiliste se produit dans le processus de racine unitaire. Même s'il existe un taux de dérive défini, il est impossible de le distinguer.

Distribution de y avec racines unitaires le dernier jour

Dans le processus de racine unitaire, l'écart type augmente à mesure que le nombre d'échantillons augmente. La valeur moyenne de y avec une taille d'échantillon de n est zéro à c = 0, mais la distribution augmente à mesure que l'écart-type $ \ sqrt {n \ cdot \ sigma ^ 2} $ augmente.

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

La distribution est différente de la précédente. Plus la taille de l'échantillon est grande, plus la distribution est large. Essayons de tracer une série.

plt.plot(s)

image.png

Dans ce cas, le niveau de chaque échantillon change par étapes. Ce changement est dû à un choc stochastique $ u_t $. Ce à quoi il a été ajouté est tracé. Une série de dérives stochastiques est tracée. Cela peut former une grande tendance. Il correspond à la zone aux deux extrémités de la carte de répartition de n = 100. Quelle plage dépend de la définition de la tendance. Cette tendance est une tendance probabiliste. Cela se distingue clairement des tendances déterministes créées par la stabilité des tendances.

Processus de retour automatique de premier ordre

Le processus de retour automatique de premier ordre peut être réécrit comme suit. 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) Ici $ y_i = c + d \ cdot i $

Lorsque $ \ rho = 1 $, ce sera un processus de marche aléatoire.

#Initialisation
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)Génération de processus
#x0 est le prix initial
#nsample est le nombre de périodes dans un essai
# rho=1: random walk
## c=c1-rho est la dérive, c1=Lorsqu'il est 1, il n'y a pas de dérive Marche aléatoire d est une tendance définie
# rho<Modèle à retour automatique 1 au 1er ordre
##Les effets de x0, c et d disparaissent avec le temps
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

Distribution de y le dernier jour du processus AR (1)

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 Dans le modèle de retour automatique du premier ordre, vous pouvez voir que lorsque $ \ rho $ devient inférieur à 1, il erre autour de la valeur initiale. A mesure que $ \ rho $ s'approche de 1, il devient difficile de le distinguer d'une marche aléatoire. La figure suivante est pour n = 7. Dans environ 7 jours, la propagation de tout $ \ rho $ sera similaire.

image.png

Détection du pouvoir (Introduction aux tests d'hypothèses statistiques par des modèles de statistiques)

Dans le test d'hypothèse, l'hypothèse nulle et l'hypothèse alternative sont faites. L'hypothèse nulle consiste en une population. Comme la population est généralement inconnue, une estimation peut être utilisée. Les spécimens sont traités comme extraits de cette population. Dans la figure suivante, la répartition de la population est indiquée par la ligne pointillée bleue. L'échantillon donne la distribution d'échantillon de cette statistique. Il est indiqué par une ligne orange. Comparez les deux. Si les deux distributions sont proches, il est difficile de juger que les deux distributions sont différentes, et si elles sont éloignées, on peut juger qu'elles sont différentes. image.png Ce jugement peut être fait en utilisant le niveau de signification. Le niveau de signification $ \ alpha $ est la zone la plus à droite de la population. 10%, 5%, etc. sont utilisés. La zone à droite de la ligne verticale bleue dans la distribution bleue. La ligne verticale bleue montre la limite où l'hypothèse nulle est rejetée. La zone à l'intérieur de cette ligne est appelée zone d'acceptation et la zone à l'extérieur de cette ligne est appelée zone de rejet.

Une fois l'échantillon obtenu, des statistiques telles que la moyenne et la variance peuvent être calculées. La valeur p est calculée à partir de la distribution de la population, qui est la probabilité d'être rejeté lorsque l'hypothèse nulle est rejetée avec la statistique comme limite inférieure. Ensuite, sur la figure, il s'agit de la zone à droite de la ligne verticale rouge. Si cette valeur est inférieure au critère de signification, on estime que l'échantillon (statistique) a été obtenu à partir d'une distribution différente de la population (population). Sur la figure, la ligne verticale rouge se trouve à gauche de la ligne verticale bleue, mais si la valeur moyenne de la distribution de l'échantillon orange est à droite de la ligne verticale bleue, il est clair que sa valeur p est inférieure au niveau de signification. image.png La capacité de rejeter correctement une fausse hypothèse nulle est appelée détection. La puissance de détection $ \ beta $ correspond à la zone à droite de la distribution orange des lignes verticales bleues. Si la partie $ \ beta $ est grande, la probabilité que l'hypothèse alternative soit correcte augmentera lorsque l'hypothèse nulle est rejetée. image.png Afin d'augmenter ce $ \ beta $, il est possible que les deux distributions soient éloignées, ou que les positions des centres des distributions soient significativement différentes. L'autre est la grande taille du spécimen. Ensuite, la forme de la distribution sera rétrécie, et si la taille de l'échantillon est suffisamment grande pour considérer la distribution comme une ligne, alors les deux lignes droites seront comparées. image.png

Test de Dicky-Fuller

Un processus AR (1) simple $ y_t= \rho y_{t-1}+u_t$ Et soustrayez $ y_ {t-1} $ des deux côtés y_t-y_{t-1}=(\rho-1)y_{t-1}+u_t ça ira. Comme $ \ delta = \ rho-1 $, l'hypothèse nulle est qu'il y a une racine unitaire dans le modèle AR. H_0:\delta=0 H_1: \delta<0 Ce sera.

Il existe trois types.

  1. Marche aléatoire sans dérive \Delta y_t=\delta y_{t-1}+u_t
  2. Marche aléatoire avec dérive \Delta y_t=c+\delta y_{t-1}+u_t
  3. Marche aléatoire avec dérive + tendance temporelle \Delta y_t=c+c\cdot t+\delta y_{t-1}+u_t Ce sont des différences de premier ordre sur le côté gauche, ce qui signifie que les différences ont des dérives ou des tendances temporelles.
# Dicky-Test plus complet
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)

Test de marche de Landau

Marche aléatoire sans dérive

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

Marche aléatoire avec dérive

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

Marche aléatoire + dérive + tendance déterministe

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

Si la marche aléatoire doit être rejetée

Modèle AR (1) sans dérive

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

Modèle AR (1) avec dérive

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

Modèle AR (1) + dérive + tendance déterministe

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

Nous pouvons voir que nous devons utiliser un modèle avec les spécifications correctes, que la taille de l'échantillon doit être grande et que ρ doit être éloigné de 1 lors du rejet de l'hypothèse nulle.

À propos de l'influence du bruit

Utilisons la distribution de Cauchy pour générer un gros bruit de queue et utilisons-la pour voir comment elle affecte la série chronologique. Regardons d'abord la distribution de la distribution de Cauchy.

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

Vous pouvez voir que la grosse queue est forte. Considérez cela comme du bruit.

Ensuite, formez la somme des nombres aléatoires de la distribution de Koshi, et regardez la distribution de la dernière valeur.

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 Après tout, vous pouvez voir que l'écart est plus grand que la distribution normale.

Ensuite, créons une fonction AR1 qui utilise la distribution de Cauchy.

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

Le test DF est effectué à l'aide de la série temporelle d'auto-retour du premier ordre créée.

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

En raison de l'influence de la grosse queue, il s'agit dans ce cas d'une tendance temporelle avec dérive, mais lorsque la taille de l'échantillon est de 250, l'hypothèse nulle ne peut être rejetée au niveau de signification de 5%. Lorsque rho = 0,85, le résultat est le suivant.

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

Lorsque la taille de l'échantillon est de 250, l'hypothèse nulle peut être rejetée au niveau de signification de 5%.

Ensuite, coupons la partie grosse queue de la distribution de Cauchy.

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

Le résultat était différent du précédent. Seules les tendances temporelles avec dérive peuvent rejeter l'hypothèse nulle, mais le niveau de signification est de 10%. Lorsque 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

Le résultat est.

Ensuite, essayez le test Dicky Fuller étendu.

\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 En augmentant la différence de premier ordre du retard, les effets structurels tels que l'autocorrélation sont supprimés.

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

Il ne peut pas être jugé correctement avec rho = 0,90. Si rho = 0,85, le résultat sera

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

J'ai pu rejeter la tendance dérive + temps.

Détection de la dérive et de la tendance temporelle lorsque l'écart type est nul ou faible

#dérive uniquement, l'écart type est nul
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

L'hypothèse nulle est rejetée du tout lorsqu'il n'y a qu'une dérive.

#dérive et tendance temporelle, l'écart type est nul
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) + dérive + tendance temporelle rejette l'hypothèse nulle quelle que soit la taille de l'échantillon

Ensuite, ajoutez un écart type à la dérive.

# 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

Ensuite, ajoutez un écart type à la tendance dérive + temps.

# 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

référence

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

Racine de l'unité

[Test Dicky-Fuller](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

Test de Dicky Fuller par des modèles statistiques
Jugement des nombres premiers par Python