R-TimeSeries

1 Comandos R para análises de séries temporais

1.1 Introdução

O R é um software livre de larga utilização em Estatística. Particularmente para Séries temporais, costuma-se, no entanto, utilizar o S-plus como software principal de análise. Isto se dá por conta de sua facilidade de utilização, poder de ajuste dos modelos e, principalmente, pela larga quantidade de modelos implementados nesta plataforma. Talvez a principal ferramenta do S-plus neste campo seja o pacote S+FinMetrics, que possui um grande número de modelos e funções de aplicações em series temporais financeiras. O Objetivo deste texto é, de forma simples, descrever de que forma alternativa podemos analisar séries temporais utilizando o software livre R. Como o R é de código aberto ele tambem possui um amplo número de implementações que nos auxiliam em tal tarefa. Portanto, diferente do S-plus - em que carregamos apenas a biblioteca FinMetrics para ter acesso a quase que a totalidade das funções S implementadas para séries temporais - no R precisaremos carregar diversos pacotes diferentes, e utilizar aquele que for mais conveniente para ajustar os modelos desejados. Outro ponto importante é que nem sempre a solução e única, no sentido de que possuímos a mesma função implementada (não necessariamente da mesma forma), em pacotes diferentes. É o caso do GARCH, por exemplo, que podemos ajustar utilizando o pacote tseries,rugarch ou ainda o fGarch.

1.2 Pacotes (bibliotecas, packages, ou libraries) utilizados neste texto:

library(stats)    #  Modelos ARIMA
library(e1071)    #  Assimetria e Curtose
library(fGarch)   #  Modelos de heterocedasticidade
library(tseries)  #  Testes de RU  
library(rugarch)  #  Modelos EGARCH univariados
library(rmgarch)  #  Modelos GARCH multivariados
library(MTS)      #  Modelos VARMA
library(FitAR)    #  Ajuste de modelos AR. Útil para plotar Estatísticas de Ljung-Box vs Lags
library(evir)     #  Calculo de Valores Extremos para VaR
library(urca)     #  Testes de Co-integração e de raízes unitárias

# Outras bibliotecas interessantes:

library(astsa)          # Shumway and Stoffer, 2011
library(finTS)          # Tsay, 2005
library(arfima)         # Fractional ARIMA Modeling
library(GECStableGarch) # ARMA-GARCH/APARCH models with GEV and stable distributions
library(gogarch)        # Generalized Orthogonal GARCH - GO-GARCH
library(gss)            # General smoothing Splines
library(evd)            # Calculo de Valores Extremos para VaR
library(fExtremes)      # Calculo de Valores Extremos para VaR
library(ismev)          # Calculo de Valores Extremos para VaR

library(ccgarch)
library(fArma)
library(fBasics)
library(fExtremes)
library(fUnitRoots)
library(tsDyn)
livrary(vars)

1.3 Análise exploratória e manipulação:

Para transformar o argumento principal (data) em um objeto da classe ts:

ts(data = NA, start = 1, end = numeric(), frequency = 1,
   deltat = 1, ts.eps = getOption("ts.eps"), class = , names = )

Para calcular estatísticas descritivas da série:

summary(object, ..., digits = max(3, getOption("digits")-3))

Desvio padrão, Assimetria e Curtose:

 sd(x, na.rm = FALSE) #Desvio Padrão
 skewness(x, na.rm = FALSE, type = 3) # Assimetria
 kurtosis(x, na.rm = FALSE, type = 3) # Curtose
 # OBS: Desvio Padrao. Se a série possuir valores missing, coloque na.rm=TRUE

Variância, Covariância e Correlação

var(x, y = NULL, na.rm = FALSE, use)

cov(x, y = NULL, use = "everything",
    method = c("pearson", "kendall", "spearman"))

cor(x, y = NULL, use = "everything",
    method = c("pearson", "kendall", "spearman"))

Para tomar a diferença da série, i.e., cria um novo objeto \(Y_t=X_t-X_{t-1}\):

diff(x, lag = 1, differences = 1, ...)

Quantis observados de X:

quantile(x, ...)

1.4 Gráficos

Para imprimir gráficos em um painel com a linhas e b colunas, com a e b inteiros positivos:

par(mfrow=c(a,b))

Gráficos da série:

ts.plot(..., gpars = list())  # Gráfico da Serie
hist(x, ...) # Histograma da Série

Para imprimir o gráfico da distribuição Normal sobre a figura do histograma, faça:

curve(dnorm(x, mean = mean(x), sd=sd(x)),add=T)

O mesmo vale se desejar imprimir a curva de outras distribuições, substibuindo dnorm pela distribuição desejadas e modificando os parâmetros da distribuição de forma conveniente.

Gráfico Quantil-Quantil para distribuição Normal:

qqnorm(y, ylim, main = "Normal Q-Q Plot",
       xlab = "Theoretical Quantiles", ylab = "Sample Quantiles",
       plot.it = TRUE, datax = FALSE, ...) 

Comando para adicionar linha central no QQplot:

qqline(y, datax = FALSE, distribution = qnorm,
       probs = c(0.25, 0.75), qtype = 7, ...) 

Esta função imprime graficos QQ não necessariamente Normais. Neste exemplo, temos um grafico QQ para um distribuição t com 5 graus de liberdade:

qqplot(y, rt(300, df = 5))

Exemplos de modificações nos gráficos de ACF e PACF:

acf0<-function(x,na.action=na.pass,lags=1:36,...){
  #lags:horizonte dos lags
  plot(acf(x,na.action=na.action,plot=F)[lags],lwd=5,lend=2,col="darkblue",...)
  axis(1, tck = 1, col = "lightgrey", lty = "dotted")
}

name="Série de log-Retornos X"
par(mfrow=c(1,2))
acf0(x, na.action = na.pass,main=paste("ACF:",name))
pacf(x, na.action = na.pass,main=paste("PACF:",name),lwd=5,lend=2,col="darkblue")

Exemplo de saída dos gráficos:

1.5 Testes:

1.5.1 Teste t:

t.test(x, y = NULL,
       alternative = c("two.sided", "less", "greater"),
       mu = 0, paired = FALSE, var.equal = FALSE,
       conf.level = 0.95, ...)

1.5.1.1 Box-Pierce | Ljung-Box (pacote stats):

Box.test(x, lag = 1, type = c("Box-Pierce", "Ljung-Box"), fitdf = 0)

1.5.1.2 Box-Pierce | Ljung-Box (pacote FitAR):

#LBQPlot(res, lag.max = 30)
### res é a série que temos interesse de realizar os testes
### lag.max é o número de lags que se deseja imprimir o gráfico
library(FitAR)
### Exemplo com 1000 observações geradas de uma normal(0,1)
res<-rnorm(1000, mean=0,sd=1)
LBQPlot(res,15)

1.5.2 Testes de Normalidade:

#### Utilizando a library `tseries`:
jarque.bera.test(x)

#### Utilizando a library `fBasics`:
ksnormTest(x, title = NULL, description = NULL)
jbTest(x, title = NULL, description = NULL)
shapiroTest(x, title = NULL, description = NULL)
normalTest(x, method = c("sw", "jb"), na.rm = FALSE)
jarqueberaTest(x, title = NULL, description = NULL)
dagoTest(x, title = NULL, description = NULL)
adTest(x, title = NULL, description = NULL)
cvmTest(x, title = NULL, description = NULL)
lillieTest(x, title = NULL, description = NULL)
pchiTest(x, title = NULL, description = NULL)
sfTest(x, title = NULL, description = NULL)

1.5.3 Testes de Raizes Unitárias

1.5.3.1 Testes de Dickey-Fuller utilizando o pacote tseries

### Testes de Dickey-Fuller
adf.test(x, lags = 1, type = c("nc", "c", "ct"), title = NULL, 
    description = NULL)
# 'lags' define o número de raízes unitarias a serem testadas
# 'type' Define o tipo de teste a ser realizado:
# 'nc' : Teste de Dickey-Fuller usual (Estatistica Tau)
# 'c'  : Teste de Dickey-Fuller Aumentado, com constante (Estatistica Tau_mu)
# 'ct' : Teste de Dickey-Fuller Aumentado, com constante e tendencia (Estatistica Tau_tau)

1.5.4 Teste de Philips-Perron (usando o pacote tseries)

pp.test(x, alternative = c("stationary", "explosive"),
        type = c("Z(alpha)", "Z(t_alpha)"), lshort = TRUE)
# 'alternative' define a Hipótese alternativa do teste
# 'stationary' : define H1: Phi < 1
# 'explosive'  : define H1: Phi > 1
#

1.6 Ajuste de modelos Univariados

1.6.1 Modelo ARIMA(p,d,q) com o pacote stats

arima(x, order = c(p,d,q))

Alternativamente, é possível utilizar as funções arma do pacote tseries, as funções do pacote astsa ou do pacote tseries. É grande a quantidade de pacotes que tem os algoritmos para modelos ARIMA implementados.

1.6.2 Modelo \(SARIMA(p,d,q)\times(P,D,Q)_s\)

arima(x, order = c(p,d,q),
      seasonal = list(order = c(P,D,Q), period = s))

1.6.3 Modelo ARIMA(p,d,q) com o pacote fArma

O pacote fArma possui a função ArmaFit, que também possibilita o ajuste de modelos ARMA incompletos. A sintaxe padrão é do tipo:

armaFit(x~arma(p,q), data, method = c("mle", "ols"), include.mean = TRUE, fixed = NULL, title = NULL, description = NULL, ...)

• (p,q) é a ordem do modelo;
  
• include.mean = TRUE adiciona a constante no modelo;
  
• Em method, se escolhe o método de estimação via Max. Verossimilhança ou Mínimos quadrados ordinários;
  
• fixed é um vetor que contém NA na posição dos pâmetros a serem estimados. Por exemplo, se quisermos estimar um ARMA(p,q), com constante, teremos p+q+1 coeficientes a serem estimados e, consequentemente fixed será um vetor com p+q+1 posições, nas quais as p primeiras serão referentes aos coeficientes da parte AR, as posições p+1 até p+q serão referentes à parte MA e a posição p+q+1 será referente à constante. Para ilustrar, se quisermos um ARMA(3,3) com constante, como definido abaixo e sem os coeficientes de \(\phi_2 X_{t-2}\) e \(\theta_1 \epsilon_{t-1}\). Entao o meu vetor fixed deve ser fixed=c(NA,0,NA,0,NA,NA,NA). Para ficar mais claro, pode-se declarar: fixed=c(c(NA,0,NA),c(0,NA,NA),NA).

\[ (1-B)^3(X_t-\phi_0)=(1-B)^3\epsilon_t \]

Alternativamente, pode-se utilizar a função declarada a seguir. Veja a explicação sobre a estrutura do parâmetro coefs na declaração da função.

Fit_Arima<-function(formula,data=NULL,order,
                       coefs=NULL, include.mean=TRUE,...){
#Para definir os elementos 1-a-1 que serão estimados no modelo, defina:
#O parametro coefs é uma lista com dois argumentos:
# $ar : vetor contendo os coefs da parte AR que devem entrar no modelo
# $ma : vetor contendo os coefs da parte MA que devem entrar no modelo
# Por exemplo, se quisermos um modelo Y_t = a_5 Y_{t-5} + a_9 Y_{t-9} - b_2 e_{t-2} - b_4 e{t_4}
# coefs deverá ser definido por: coefs = list(ar=c(5,9),ma=c(2,4))

# Montando a estrutura do argumento 'fixed' para ser passado para a funcao 'arima'
fixed=NULL
  if(!is.null(coefs)){
    fixed<-rep(0,times=sum(order,include.mean))
    if(include.mean==TRUE){
    fixed[c(coefs$ar,order[1]+coefs$ma,length(fixed))]<-NA
    }
    else {
    fixed[c(coefs$ar,order[1]+coefs$ma)]<-NA
    }
  }
  fit<-armaFit(as.formula(formula), data=data, fixed=fixed, include.mean=include.mean,...)
  return(fit)
}

1.6.4 Outra alternativa para ajustar o modelo ARMA incompleto

inc.arima<-function(x,order=NULL,fixed=NULL,lag=NULL,name="Log-Retornos de X"){
if(!is.null(lag)){ 
  # vetor de zeros do tamanho do numero de parametros a serem estimados
    fixed<-rep(0,times=sum(order,1)) 
    # Preenchendo com 'NA' as posicoes dos parametros que desejamos que sejam estimados
    fixed[c(coefs$ar,order[1]+coefs$ma,length(fixed))]<-NA 
      fit<-arima(x,order=order,fixed=fixed) #Ajustando o modelo
      return(fit)
}

# Se o parametro `ccoefs` não for passado, entao ajusto o modelo completo.
fit<-arima(x,order=order)
return(fit)
}

1.6.5 Ajuste de modelos ARMA(p,q)-GARCH(m,k) utilizando o pacote fGarch

fit<-garchFit(x ~ arma(p,q)+garch(m, k), data =  banco_de_dados, 
               cond.dist = c("std"),include.mean=TRUE)
               
# O parametro 'cond.dist' define a distribuição dos erros a ser utilizada: 
# 'norm'  : Erros normais
# 'std'   : Erros t-student
# 'ged'   : Dists de Erros Generalizados

# Para imprimir detalhes do ajuste, como tabela de analise de variancia 
# e testes de Ljung-Box, aplique 'summary' ao modelo ajustado:
summary(fit)          

1.6.6 Modelo EGARCH(1,1) utilizando o pacote rugarch

# Neste pacote, é usual especificar todos os parametros inicialmente e informá-los 
# todos de uma vez à funcão de ajuste:
spec=ugarchspec(variance.model=list(model="eGARCH", garchOrder=c(1,1)),
mean.model=list(armaOrder=c(1,1), arfima=FALSE),
distribution.model="std")       # Especificacao dos parâmetros
fit=ugarchfit(data=x,spec=spec) # Ajuste do modelo
summary(fit)                    # Resultados detalhados do ajuste

1.6.7 Modelos TGARCH(m,k) utilizando o pacote fGarch

O modelo TGARCH é um caso específico do modelo APARCH considerando delta=1 (Laurent,2004). Utilizamos o pacote fGarch para ajustar um modelo APARCH(1,1),\(\delta\)=1, que é equivalente a um TGARCH(1,1):

require(fGarch) 
fit_TGARCH<-garchFit(formula = ~ arma(p,q)+aparch(1, 1), data =  x,
                   cond.dist = c("std"),include.mean=TRUE,
                   include.delta=TRUE,delta=1)

1.6.8 Análise de resíduos para o modelo GARCH

Para avaliar se o modelo ARMA-GARCH ajustado está captando as dependências anteriores do nível e da volatilidade, usualmente se constrói os gráficos ACF e PACF dos resíduos e quadrados dos resíduos. Ainda, alternativamente aos pacotes fGarch e rugarch, o pacote tseries também possui uma função garch. Aqui, a única ressalva é com relação às funções do pacote fGarch, que produzem resíduos não-padronizados, diferente dos outros dois pacotes. Logo, suponha que temos um ajuste de um modelo via o pacote fGarch, digamos fit. Para obter resultados equivalentes, comparáveis e interpretáveis, a sintaxe correta para se obter os resíduos de um modelo ajustado em fit é at=residuals(fit_norm_cemig_at,standardize = T)).
A partir, então de at, pode se construir as análises de \(a_t\) e \(a_t^2\) da forma usual, isto é acf(at),pacf(at),acf(at^2),pacf(at^2).

1.6.9 Previsões

Tanto para os modelos ARMA quanto para os modelos da família GARCH apresentados aqui, a sintaxe para se obter previsões é:

predict(fit,n.ahead=h)

em que h é o horizonte de previsão e fit é um objeto contendo o modelo ajustado.

1.7 VaR: Value-at-risk

1.7.1 Metodologia RiskMetrics

A seguinte função pode ser facilmente programada no R para o cálculo do VaR:

RMVaR<-function(rt, lambda = 0.94, k = 1, long = TRUE, alpha = 0.95,range.optim=seq(0.9,1,0.001),plot=FALSE){
  #### rt é a série de log-retornos
  #### se long == TRUE, entao estamos no caso de uma posicao comprada. Se long == FALSE, entao a posicao será considerada vendida
  #### lambda é o parâmetro de suavização
  #### 'alpha' é a probabilidade fixada da distribuição Normal padrao - Lembrete: o RiskMetrics supõe erros normais.
  #### k é o horizonte para o qual se deseja calcular o VaR
  #### Se lambda == NULL, entao o programa executa a rotina para todos os valores em range.optim=seq(0.9,1,0.001)
  #### Se lambda == NULL e plot == TRUE, entao o programa imprime a soma de quadrado dos resíduos para os lambda em range.optim
 
  # Se long=T, entao estou em uma posicao comprada e o VaR é definido na cauda inferior, i.e., lower.tail=T.
  # Se long=F, entao o VaR é calculado na cauda superior e lower.tail = F.
  
  Z=qnorm(0.05,lower.tail=long)
  r2t<-rt^2
  T=length(r2t)
  SQR<-numeric()
  VaR<-numeric()
  sigma2hat<-numeric()
  sigma2hat[1]<-r2t[1]
  
  if(!is.null(lambda)){ # Entao nao ira para a otimizacao de lambda
    for(t in 2:T){
      sigma2hat[t]<-sigma2hat[t-1]*lambda+(1-lambda)*r2t[t-1]
      #print(paste("hello 24-",t))
    }
    lambda_ = lambda
    S2t1<-sigma2hat[T]*lambda+(1-lambda)*r2t[T]
    VaR<-sqrt(S2t1*k)*Z
    SQR<-sum((sigma2hat-r2t)^2,na.rm = TRUE)
  }
  else {
    if(!is.null(range.optim)){
      for( i in 1: length(range.optim)){
        tmp<-RMVaR(rt, lambda = range.optim[i], k = k, long = long, alpha =alpha)
        VaR[i]<-tmp[[1]]
        SQR[i]<-tmp[[2]]
      }
      if(plot==TRUE){
        plot(range.optim,SQR, type="l", col="darkblue",lwd=2, main="RiskMetrics\n Sum of Squared Residuals",xlab="lambda",ylab="SQR")
        abline(h=min(SQR),v=range.optim[SQR==min(SQR)],lty=2,col="gray")
      }      
    }
    else {
      return("Error! You must specify 'lambda' or 'range.optim' ")
    }
  lambda=range.optim[SQR==min(SQR)]
  lambda_=range.optim
  print(paste("RiskMetrics - The Minimum SQR is: ",round(min(SQR),4)," at lambda =",lambda))   
  }
  return(list(VaR=VaR,SQR=SQR,lambda=lambda_))
}

Exemplo de Utilização utilizando dados e log-retornos do Ibovespa:

load("VaR.RData") 
RMVaR(r_ibov,long=F,k=30) #Cálculo do VaR para periodo de 30 dias com o lambda padrão, sugerido pelo RiskMetrics, i.e., lambda=0.94

## $VaR
## [1] 0.1962522
## 
## $SQR
## [1] 0.01421295
## 
## $lambda
## [1] 0.94

RMVaR(r_ibov,long=F,k=30,lambda=0.88) #Cálculo do VaR para periodo de 30 dias com lambda = 0.88

## $VaR
## [1] 0.2004856
## 
## $SQR
## [1] 0.01417891
## 
## $lambda
## [1] 0.88

Optim<-RMVaR(r_ibov,long=F,k=30,lambda=NULL,range.optim=seq(0.8,1,0.001),plot=T) # Cálculo do VaR para periodo de 30 dias com lambda selecionado iterativamente de acordo com a menor soma de quadrados dos residuos

## [1] "RiskMetrics - The Minimum SQR is:  0.0142  at lambda = 0.901"

1.7.2 Função para imprimir gráfico de quantis empíricos

qtEmpirico<-function(x,p,long=TRUE,plot=T,nbreaks=30,print=TRUE,title="Empirical VaR"){
  dnst<-density(x)
  qt<-quantile(x,p*long+(1-p)*!long)
  if(plot==T){
    hist(x,breaks=nbreaks,col="lightgrey",border="white",freq=F,
         ylim=c(0,max(dnst$y)),xlim=c(min(x),max(x)),main=title)
    par(new=T)
    if(long != TRUE){ ### Posicao Vendida
      plot(dnst$x[dnst$x>qt],
           dnst$y[dnst$x>qt],type="h",
           ylim=c(0,max(dnst$y)),xlim=c(min(x),max(x)),
           xaxt="n",yaxt="n",lwd=1,col="blue",ylab="",xlab="",main="")
    }
    else{ ### Posicao comprada
      plot(dnst$x[dnst$x<qt],
           dnst$y[dnst$x<qt],type="h",
           ylim=c(0,max(dnst$y)),xlim=c(min(x),max(x)),
           xaxt="n",yaxt="n",lwd=0,col="blue",ylab="",xlab="",main="")  
    }
    par(new=T)
    plot(dnst,xaxt="n",yaxt="n",
         ylim=c(0,max(dnst$y)),xlim=c(min(x),max(x)),
         lwd=1,col="darkblue",ylab="",xlab="",main="")
  }
  if(print==TRUE){
  print(paste0("Probability: ",p))
  print(paste0("Position long?: ",long))
  print(paste0("Empirical Quantile: ",round(qt,5)))
  }
  return(qt)
}
qtEmpirico(rnorm(500),0.05,long=T)

## [1] "Probability: 0.05"
## [1] "Position long?: TRUE"
## [1] "Empirical Quantile: -1.75689"

##        5% 
## -1.756889

1.7.3 Distribuição de Valores Extremos

O Pacote evir é análogo do R ao pacote evis do S-plus Assim, para estimar os parâmetros da distribuição de Valores extremos:

\[ \begin{align} VaR & = \mu_n - \frac{\sigma_n}{\xi_n}\left\{1-\left(-n log(1-p)\right)^{\xi_n}\right\} , \xi_n\neq 0 \\ & = \mu_n - \sigma_n log\left(-n log(1-p)\right) , \xi_n = 0 \\ VaR(k) & = k^{\xi}VaR \end{align} \]

utilizamos o comando gev:

gev(data,block)

• data é o objeto contendo os dados que se deseja estimar a distribuição;
  
• block é o tamanho da janela utilizada para se calcular os máximos locais.

Pode-se ainda utilizar de forma adaptada a função qtEmpirico(.) apresentada anteriormente para se calcular o VaR via a distribuição de valores extremos estimada aqui.

1.8 Modelos Multivariados

Para toda esta seção, considere x um data.frame cotendo as séries de interesse, todas elas com mesmo tamanho;

1.8.1 Função de correlação cruzada

Cálculo das matrizes de covariancia e representação pictórica utilizando o pacote MTS:

ccm(x,level=T)

• Se, level=F, então ccm imprime somente a representação pictórica, sem os valores das matrizes.

1.8.2 Testes para escolha de modelos

Os comandos VARorder e VMAorder calculam as estatísticas AIC, BIC e HQ para diversas ordens dos modelos VAR(p) e VMA(q) respectivamente.

VARorder(x)

VMAorder(x)

1.8.3 Ajuste dos Modelos

1.8.3.1 VAR(p) - Modelo AR multivariado com ordem p

VAR(x,p=p,include.mean=T)

1.8.3.2 VMA(q) - Modelo MA multivariado com ordem q

VMA(x,p=q,include.mean=T)

1.8.3.3 VAR(p) incompleto - Modelo AR multivariado com ordem p incompleto

VARs(x,p=p,include.mean=T,lags=coefs)

1.8.3.4 VMAs(q) incompleto - Modelo MA multivariado com ordem q incompleto

VMAs(x,p=q,include.mean=T,lags=coefs)

1.8.3.5 VARMA(,p,q)

VARMA(x,p=p,q=q,include.mean=T)

1.8.3.6 Qualidade de Ajuste dos modelos

MTSdiag(fit,adj=12)

• Fit é um objeto contendo o ajuste do modelo VAR, VMAR ou VARMA de interesse. MTS calculará as estatísticas de Ljung-Box e Portmanteau para os resíduos do modelo. Para detalhes, veja Tsay(2014).

1.8.3.7 Detalhes:

• Em todos os 3 modelos acima, include.mean=T adiciona a constante no modelo. Se include.mean=F, então o modelo estimado não considerará a constante;
• fixed é um vetor com a mesma estrutura daquela apresentada para os modelos do pacote fArma, isto é, ele deve ser do tamanho do número de parametros a serem estimados e deve conter NA nas posições dos coeficientes de interesse e zero nos que se deseja excluir.
• Para os comandos VARs e VMAs, é possivel omitir alguns coeficientes, como tinhamos no caso dos modelos ARMA incompletos. Neste caso, coefs é um vetor com as posições que se tem interesse que sejam ajustadas, e.g., se queremos um VAR(3) sem o coeficiente de \(\underline{X}_{t-2}\), executamos: VARs(x, p=3,include.mean=TRUE,lags=c(1,3))

1.9 Co-integração

Os testes de co-integração desta seção serão realizados com o pacote urca do R.

1.9.1 Teste de Phillips & Ouliaris:

ca.po(X, demean = c("none", "constant", "trend"),lag = c("short", "long"))

X é um objeto do tipo data.frame contendo as duas séries as quais se tem interesse de testar a existência de cointegração.
Defina demean =constant para considerar as séries do modelo com constante (drift). Em outras palavras, utilize esta opção se as séries envolvidas possuírem média constante diferente de zero. A opção demean =trend considera o caso com séries com tendência e drift e o caso demean =none contempla a situação de séries ambas sem tendência e drift.
lag é a opção que define a ordem máxima do modelo VAR ajustado. Se lag='short' então a ordem máxima do VAR ajustado será \(\left[4\left(\frac{T}{100}\right)^{0.25}\right]\). Se lag='long' então teremos a ordem máxima igual a \(\left[12\left(\frac{T}{100}\right)^{0.25}\right]\), em que T indica o comprimento a série.
A forma padrão da função ca.po utiliza a estatística \(P_u\) (Veja Philips e Ouliaris, 1990) para realizar o teste. Desta forma, o teste rejeita H: As séries não são co-integradas para valores grandes da estatistica Pu.
Exemplo(Morettin, 2013):

\[ \begin{alignat}{1} X_{1t} & = & \beta_2 X_{2t} + \beta_3 X_{3t} + u_t \\ X_{2t} & = & X_{2,t-1} + v_t \\ X_{3t} & = & X_{3,t-1} +w_t \\ \underset{\sim}{\beta}^T & = & (1,-\beta_2,-\beta_3)=(1,-0.5,-0.5) \\ u_t & = & 0.75 u_{t-1}+\epsilon_t \nonumber \end{alignat} \]

# Simulando o sistema
beta=c(1,-0.5,-0.5)

u<-arima.sim(n = 250, list(ar = c(0.75),ma=NULL), sd = 0.5)
v<-rnorm(n=250, mean=0, sd=0.5)
w<-rnorm(n=250, mean=0, sd=0.5)

X2<-numeric()
X3<-numeric()
for(t in 1:250){
  if (t==1) {
    X2[t] = 0 + v[t] # Inicializando com X2[0]=0
    X3[t] = 0 + w[t] # Inicializando com X3[0]=0
  }
  else {
    X2[t] = X2[t-1] + v[t]
    X3[t] = X3[t-1] + w[t]
  }
}

X1= -beta[2]*X2-beta[3]*X3+u # Relação de cointegração
X<-cbind(X1,X2,X3) # Cria a matriz de dados

library(urca)
cointest_X<-ca.po(X, demean="constant") # Teste de P&O para o caso com constante
summary(cointest_X)

## 
## ######################################## 
## # Phillips and Ouliaris Unit Root Test # 
## ######################################## 
## 
## Test of type Pu 
## detrending of series with constant only 
## 
## 
## Call:
## lm(formula = z[, 1] ~ z[, -1])
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.03547 -0.48997  0.01617  0.37825  1.73713 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.17036    0.05546   3.072  0.00237 ** 
## z[, -1]X2    0.40363    0.02751  14.673  < 2e-16 ***
## z[, -1]X3    0.46881    0.00958  48.935  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.6677 on 247 degrees of freedom
## Multiple R-squared:  0.9187, Adjusted R-squared:  0.9181 
## F-statistic:  1396 on 2 and 247 DF,  p-value: < 2.2e-16
## 
## 
## Value of test-statistic is: 103.2199 
## 
## Critical values of Pu are:
##                   10pct    5pct    1pct
## critical values 33.6955 40.5252 53.8731

O valor da estatística é 103.2199274 > 53.8731. Logo, rejeitamos H: Não há co-integração ao nível 1%. O vetor de co-integração é: (1, -0.4036312, -0.4688126) e a estimativa da constante é 0.1703546.

É necessário ainda avaliar os resíduos do modelo de correção de erros. O gráfico dos resíduos pode ser obtido por:

plot(cointest_X)

Se houver necessidade de se modelar estes resíduos, ou de se realizar análises adicionais, eles podem ser obtidos através da propriedade res, e armazenados em um outro objeto. Por exemplo:

resid_coint<-cointest_X@res # Armazenando os resíduos do modelo ajustado em 'resid_coint' 

Estatísticas descritivas para os resíduos do MCE ajustado no exemplo

	resid_coint
	Min. :-2.03547
	1st Qu.:-0.48997
	Median : 0.01617
	Mean : 0.00000
	3rd Qu.: 0.37825
	Max. : 1.73713

1.9.2 Teste de Johansen

Considerando o sistema multivariado dado por: \[ \underset{\sim}{X}_t=\Pi_1 \underset{\sim}{X}_{t-1}+\dots+ \Pi_k \underset{\sim}{X}_{t-k}+\mu+\Phi D_t +\epsilon_t , (t=1,\dots,T) \] O modelo vetorial de correção de erros é especificado por:
\[ \Delta\underset{\sim}{X}_t=\Gamma_1 \Delta\underset{\sim}{X}_{t-1}+\dots+ \Gamma_{k-1} \Delta\underset{\sim}{X}_{t-k+1}+\Gamma_{k}\underset{\sim}{X}_{t-k}+\Phi D_t +\epsilon_t , (t=1,\dots,T)\\ \Gamma_i=-(\underset{\sim}{I}-\Pi_1-\dots-\Pi_i),(i=1,\dots, k) \]
Em que os vetores \(\mu,D_t,\Phi,\epsilon_t,\Pi_i\) e \(\Gamma_i,i=1,\dots,k\) são especificados como em (Johansen,1988). A sintaxe para este teste, utilizando o pacote urca é:

joh_X<-ca.jo(X, type = c("eigen", "trace"), ecdet = c("none", "const", "trend"), K = 2)
summary(joh_X)

Novamente, X é do tipo data.frame. A opção type permite que se escolha qual estatística será utilizada. Utilize type="eigen" para a estatística do maior autovalor e type="trace" para a estatística do traço. A opção ecdet contempla os diversos casos possíveis de adição de constantes e tendências ao VECM. Veja (Zivot e Wang, 2006) para detalhes. K é a ordem do modelo VAR(k) a ser utilizada(Veja equação do sistema proposto acima) e é suposta igual a 2 se suprimida.
O comando summary permite visualizar os detalhes do ajuste, como a martriz de cargas, o(s) vetor(es) co-integrados e as estatísticas de teste. Novamente, rejeitamos a hipótese \(H:r\le r_0\) para valores grandes da estatística, seja ela a do traço ou a do máximo autovalor.

1.10 Referências

• Engle, R.F. and Granger, C.W.J. , Co-integration and Error Correction: Representation, Estimation, and Testing, Econometrica, Vol. 55, No.2, pp. 251-76, 1987.
  
• Hannan, E.J. and Quinn, B.G. , The Determination of the Order of an Autoregression, Journal of the Royal Statistical Society. Series B (Methodological) V41(2):190–195, Wiley, 1979.
  
• Johansen, S., Statistical Analysis of cointegration vectors, Journal of Economic Dynamics and Control, 12, 231-254, 1988.
• Laurent S., Analytical Derivates of the APARCH Model, Computational Economics, V24, Issue 1, pp 51-57, August 2004, Kluwer Academic Publishers.
  
• Morettin P.A., Econometria Financeira, Blucher 2ed. Sao Paulo, 2013.
  
• Morettin P.A. e Toloi, C.M.C, Análise de Séries Temporais, Blucher, Sao Paulo, 2006.
  
• Phillips, P.C.B. and Ouliaris, S., Asymptotic Properties of Residual Based Tests for Cointegration, Econometrica, Vol. 58, No. 1, 165–193, 1990.
  
• Shumway R.H. and Stoffer D.S., Time Series Analysis and Its Applications: With R Examples, 3rd ed. Springer Texts in Statistics, 2011.
• Taylor, J.W., A Quantile Regression Neural Network Approach to Estimating the Conditional Density of Multiperiod Returns, Journal of Forecasting, 19:299–311, Wiley and Sons, 2000.
  
• Tsay, R.S., Analysis of Financial Time Series, 2nd ed., Wiley, 2005.
• Tsay, R.S., Multivariate Analysis of Time Series with R and Financial Applications, Wiley and Sons, 2014.
• Zivot, E., Wang, J., Modeling Financial Time Series with S-PLUS, 2nd ed., Springer, 2006.