Exemplo de estimação de VAR no "R"

A idéia nesse post é "ajudar" na estimativa de um Vetor auto-regressivo usando o programa "R"....OBS: os comandos estão escritos em vermelho.

Os dados utilizados foram coletados no sistema de gerenciador de séries temporais do Bacen (link barra lateral) e no Ipeadata (link barra lateral). Para facilitar o exercício salvei o documento em minha conta no google docs, é só clicar a baixar a tabela pro teste, o nome do arquivo é "macro".
A períodicidade é mensal, de jan de 2000 a dez 2010....as variáveis já estão estacionárias (para PIB, SELIC (dessazonalizado) e CÂMBIO aplicamos a primeira diferença), IBOVESPA (Ibovespa - variação percentual mensal) e IGP-M são estacionãrias em nível.

O primeiro passo é importar os dados para o programa.....eu sempre tenho costume de, dentro do R, ir em "Arquivo" - "Mudar dir.." e escolher a pasta para importar os dados salvos.....no meu caso, salvo sempre dentro da própria pasta do R....em "c:"-"arquivos e programas" - "R" - "R-2.11.1" (salvo meus dados sempre dentro dessa pasta)....é só indicar o caminho de onde o programa tem q importar..

O comando para importar os dados é ... - macro<-read.table("macro.txt", head=T)
No R os sinais " -< e =" são equivalentes, estou importando a tabela com o nome de e nomeando como macro
para poder trabalhar com os dados attach(macro)
Caso você queira ter algumas informações sobre as séries de dados, máx, min, média, mediana summary(macro)

Depois,,.devemos carregar a library "vars"..seja manualmente em "Pacotes" - "carregar pacotes" (na barra superior), seja via comando - library("vars") - . Caso você não tenha esse pacote vá "Pacotes" - "instalar pacotes" e instale-o.

Para testar a estacionáridade da série através de um teste ADF digitamos o seguinte comando
adf1 <- summary(ur.df(macro[, "PIB"], type = "trend", lags = 2))
dei o nome de adf1 ao teste, para testar as outras séries é só substituir o "PIB" pelo nome da série entre "aspas", o teste foi realizado considerando dois lags, para modificar basta modificar o número na equação de comando. O resultado do foi o seguinte:

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1279.06230 792.86438 1.613 0.1093
z.lag.1 -1.49650 0.21126 -7.084 9.52e-11 ***
tt 21.87234 10.41879 2.099 0.0378 *
z.diff.lag1 0.25468 0.15449 1.649 0.1018
z.diff.lag2 -0.12005 0.09793 -1.226 0.2226
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Para testar a defasagem a ser utilizada pelo VAR, através dos testes AIC, BIC, HQ, FPE usamos o seguinte comando... VARselect(macro, lag.max = 8, type = "both")
usei uma defasagem de 8 lags pois os resultados em geral não passam de 3, 4 defasagens....mas para aumentar-diminuir basta mudar o número.
saída do R.....

AIC(n) HQ(n) SC(n) FPE(n)
2 2 1 2

O resultado sugere uma defasagem de ordem 2, considerando 3 dos 4 testes.

Descoberta a defasagem vamos estimar o VAR...p1ct <- VAR(macro, p = 2, type = "both")
Como são 5 variáveis e consideramos a relação endógena entre todas, mais um intecerpto e um trend temos 5 equações....para sumarizar os resultados do VAR especificos a uma variável (uma equação) .. summary(p1ct, equation = "SELIC") , para estimar as outras 4 equações basta substituir o nome da variável entre "aspas".
O resultado para SELIC é o seguinte:

Estimation results for equation SELIC:
SELIC = CAMBIO.l1 + PIB.l1 + SELIC.l1 + IBOV.l1 + IGP.l1 + CAMBIO.l2 + PIB.l2 + SELIC.l2 + IBOV.l2 + IGP.l2 + const + trend

Estimate Std. Error t value Pr(>|t|)
CAMBIO.l1 -3.586e-03 2.733e-03 -1.312 0.192063
PIB.l1 -2.222e-06 2.426e-06 -0.916 0.361540
SELIC.l1 -7.276e-01 8.877e-02 -8.197 3.56e-13 ***
IBOV.l1 -3.121e-03 1.326e-03 -2.354 0.020261 *
IGP.l1 1.427e-02 2.034e-02 0.702 0.484314
CAMBIO.l2 -1.888e-03 2.672e-03 -0.707 0.481215
PIB.l2 -4.025e-06 2.420e-06 -1.663 0.098936 .
SELIC.l2 -2.922e-01 8.615e-02 -3.392 0.000949 ***
IBOV.l2 -2.329e-03 1.470e-03 -1.585 0.115763
IGP.l2 4.173e-02 1.894e-02 2.203 0.029547 *
const -5.507e-02 2.434e-02 -2.263 0.025509 *
trend 3.498e-04 2.704e-04 1.293 0.198419

--- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Podemos ainda plotar os resíduos.... plot(p1ct, names = "SELIC")
O PLOT DOS RESÍDUOS ESTÁ NA FIGURA ACIMA, BASTA CLICAR NELA QUE O GRÁFICO MAXIMIZA

Para os testes sobre os resíduos temos...Portmanteau Test (asymptotic)
ser11 <- serial.test(p1ct, lags.pt = 16, type = "PT.asymptotic") ser11$serial
No nosso caso o resultado é....
Portmanteau Test (asymptotic)
data: Residuals of VAR object p1ct
Chi-squared = 367.9793, df = 350, p-value = 0.2440

E para teste de normalidade, JB e Kurtose,...temos o seguinte comando...
norm1 <- normality.test(p1ct)
norm1$jb.mul

Nossos resultados....

JB-Test (multivariate)
data: Residuals of VAR object p1ct
Chi-squared = 146.5162, df = 10, p-value < squared =" 6.6308," df =" 5," value =" 0.2496" squared =" 139.8854," df =" 5,"

Para o impulso resposta o comando é irf(p1ct, impulse = "SELIC", response = c("PIB"), boot = FALSE)
para mudar o impulso e a resposta basta mudar os entre "aspas", a resposta pode ser dada a mais de uma variável, para isso basta acrescentar uma vírgula e digitar a outra variável entre "aspas".
para ver o resultado basta digitar
irf
O resultado que vai aparecer para você é...
Impulse response coefficients
$SELIC
PIB
[1,] 0.00000
[2,] -526.25830
[3,] 444.91295
[4,] -19.93247
[5,] -415.43087
[6,] 216.09450
[7,] 99.90223
[8,] -163.94160
[9,] 36.27845
[10,] 48.33234
[11,] -38.06526


OBS: A estimação aqui realizada é meramente ilustrativa (estou trabalhando com dados parecidos para um artigo de uma disciplina, e ainda, existem alguns artigos rodando VAR com essas mesmas variáveis). Contudo me responsabilizo por qualquer erro que tenha cometido nesse post e caso alguem tenha alguma contribuição ou ainda reclamação favor manifestar-se.

Espero ter contribuído