Modélisation VARMA (package DSE et VARS)

[Remi Auguste]

Bonjour,

Avant tout je tiens à préciser que je suis nouveau sur R (ainsi que sur ce forum), aussi je vous demande d'être indulgent avec moi :oops: Avant de poster cette question, j'ai bien entendu recherché si une question semblable à la mienne existait (sur internet comme sur le forum), sans succès.

Voilà mon problème, j'ai une série temporelle, constituée de vecteur 2D, située dans un fichier. J'arrive sans problème à charger ces données dans une matrice. En revanche je n'arrive pas à modéliser ma série temporelle par un modèle ARMA multivarié (ou VARMA). D'après la documentation le package DSE contient les outils nécessaire pour cela. Je n'arrive hélas pas à comprendre comment passer de ma matrice à l'objet TSdata pour l'utiliser avec la fonction estVARXls.

J'ai beau lire la documentation fournie avec le package DSE, je n'arrive pas à la comprendre, notamment en ce qui concerne l'input et l'ouput de TSdata. :cry:

Jusque là voici les quelques lignes dont je sois (à peu près) sûr:

• library("dse")
• setwd("/data/writeable")
• m1 <- t(matrix(read.table("abvectors.dat")[2:3]))

La première colonne de mon fichier abvectors.dat (disponible à l'adresse: http://www2.lifl.fr/~auguster/abvectors.dat) ne contient que l'indice du vecteur.

J'apprécierais si vous pouviez me guider sur les étapes pour réaliser cette modélisation. Merci d'avance.

PS: En semaine je suis aussi disponible sur le serveur irc.freenode.net pseudo Foloex.

[Matthieu Stigler]

Bienvenue sur le forum

Pour qu'on puisse t'aider, il est bien de montrer un code, mais il devrait être reproductible!! read.table est donc strictement illégal ;-)

Code : Tout sélectionner

library("dse")
#copions l'exemple de ?TSdata
rain  <- matrix(rnorm(86*3), 86,3)
radar <- matrix(rnorm(86*2),  86,2)
mydata <- TSdata(input=radar, output=rain)

Maintenant le modèle:

Code : Tout sélectionner

 estVARXls(mydata, max.lag=2)

Maitenant, à toi d'interpréter ce que ça veut dire, je connais aps le package et l'output est asse horrible, surtout la fonction coef() donne rien de très lisible.

Si tu ne veux pas absolument utiliser un VARMA mauis plutot un VAR (plus ou moins équivalent, bcp plus utilisé en pratique à ma connaissance) utilise le package varsm qui propose en plus tests de spécification (heterosced, autocorr, normailité), tests de causalité et impulae response functions

[Remi Auguste]

Pour qu'on puisse t'aider, il est bien de montrer un code, mais il devrait être reproductible

C'est parfaitement reproductible, je fourni même mon fichier de données que j'ai utilisé...

read.data est donc strictement illégal

Ça tombe bien, je n'utilise pas read.data ...

J'avoue que même en relisant l'exemple, je ne vois pas ce que sont input et output ? Ni comment les obtenir dans mon cas.
Je vais voir si la documentation de varsm me permet de mieux comprendre. Merci :wink:

[Matthieu Stigler]

ah pardon c'est vrai c'est reproductible... y'avait meme le fichier :-(

Oui le pkg vars est bien documenté, regarde notamment dans le journal of statistical software. Sinon il y a encore quelques fonctions dans le package tsDyn qui te permettent d'estimer des VAR, mais vars est plus complet.

[Remi Auguste]

J'ai réussi à avancer dans le bon sens (enfin j'espère).
J'arrive maintenant à mettre mes données dans un objet TSdata (sans pour autant comprendre à quoi correspond output).

J'utilise maintenant ce fichier de données: www.lifl.fr/~auguster/datas.dat

Mon code est

Code : Tout sélectionner

library("dse")
m1 <- read.table("datas.dat")[2:3]
m2 <- matrix(0, dim(m1)[1], dim(m1)[2])
data1 <- TSdata(input=m1, output=m2)
model1 <- estVARXls(data1)


Tout à l'air de bien fonctionner jusqu'à estVARXls qui me créé un model inutilisable:

Code : Tout sélectionner

> model1
neg. log likelihood= 150.7059

A(L) =
1    0   
0    1   

B(L) =
1    0   
0    1   

C(L) =
0    0   
0    0 


avec les warnings:

> model1 <- estVARXls(data1)
Warning messages:
1: In lsfit(Past, outputData(data)[(max.lag + 1):N, , drop = FALSE], :
'X' matrix was collinear
2: In residualStats(r$pred, y, sampleT, warn = warn) :
The cov. matrix is singular. Working on subspace.

J'avoue que ce genre de message ne me parle pas du tout.

Avec JMulti, j'obtient un VAR(6) sur ces données, mais j'aimerais y arriver sur R. De plus, en passant data1 à la fonction bft() (que je trouve super bourine au passage :D), il s'arrête sur au lag=6, ce qui rejoint JMulti. En revanche je n'arrive toujours pas à trouver comment on estime les paramètres...

[Matthieu Stigler]

Oui comme je le disais les résultats de la fonction sont très étranges, essayes plutôt avec paquet vars:

Code : Tout sélectionner

library(vars)
VAR6_vars<-VAR(m1, p=6)
summary(VAR6_vars)
coef(VAR6_vars)

ou tsDyn:

Code : Tout sélectionner

library(tsDyn)
VAR6_tsdyn<-lineVar(m1, lag=6)
summary(VAR6_tsdyn)
coef(VAR6_tsdyn)

[Remi Auguste]

J'ai utilisé le code suivant pour estimer automatiquement l'ordre du model, estimer les paramètres et les récupérer les paramètres.

Code : Tout sélectionner

model <- VAR(m1, ic="AIC", lag.max=20)
Acoef(model)


Bon ça ne me dit toujours pas comment faire un VARMA, mais cela va me permettre d'exploiter mes données.

Merci beaucoup pour ton aide en tout cas !