Índice de Vegetação de Diferença Normalizada (NDVI) de uma célula de grade no Alasca central. NDVI é uma medida da verdura da vegetação e está relacionada à cobertura de vegetação verde, atividade fotossintética e biomassa verde. NDVI varia entre -1 e 1.

O solo descoberto e a neve geralmente tem valores de NDVI abaixo de 0.2, enquanto que as áreas vegetadas têm valores de NDVI acima de 0.2. Esta série temporal do NDVI é do conjunto de dados GIMMS (Global Inventory, Monitoring, and Modeling Studies), em Tucker et al. 2005, que fornece estimativas do NDVI a partir de observações de satélite AVHRR (Advanced Very High Resolution Radiometer).

Referência
Tucker, C.; Pinzon, J.; Brown, M.; Slayback, D.; Pak, E.; Mahoney, R.; Vermote, E.; El Saleous, N.,
An extended AVHRR 8-km NDVI dataset compatible with MODIS and SPOT
VEGETATION NDVI data. International Journal of Remote Sensing 2005, 26, 4485-4498. 

ndvi.obs = ts(read.table("http://leg.ufpr.br/~lucambio/CE017/20212S/ndvi.txt", col.names = "ndvi"), start = c(1982,1), freq = 12)
par(mar=c(4,4,1,1))
plot(ndvi.obs, pch=16, xlab="TEMPO", ylab="NDVI", type="b", col="black", main="")
abline(h=mean(ndvi.obs, na.rm = TRUE), col="red", lwd=2)
grid()

library(astsa)
monthplot(ndvi.obs, xlab="MESES", ylab="NDVI", type="l")
grid()

Observamos que, para cada mês, nos últimos anos sempre os valores de NDVI são menores do que nos primeiros anos observados. 

par(mar=c(4,4,1,1))
plot(ndvi.obs, pch=16, xlab="TEMPO", ylab="NDVI", type="b", col="black", main="")
lines(ksmooth(time(ndvi.obs), ndvi.obs, "normal", bandwidth=1), lwd=2, col=4)
grid()

Eliminando a tendência

ndvi.obs1 = ndvi.obs - ksmooth(time(ndvi.obs), ndvi.obs, "normal", bandwidth=1)$y
par(mar=c(4,4,1,1))
plot(ndvi.obs1, pch=16, xlab="TEMPO", ylab="NDVI", type="b", col="black", main="")
abline(h=mean(ndvi.obs1), col="red", lwd=2)
grid()

Observando a estrutura de autocorrelção

lag1.plot(ndvi.obs1, 12)

acf2(ndvi.obs1)

##      [,1]  [,2]  [,3]  [,4]  [,5]  [,6]  [,7]  [,8]  [,9] [,10] [,11] [,12]
## ACF  0.64  0.35 -0.01 -0.37 -0.67 -0.74 -0.65 -0.34 -0.01  0.36  0.62  0.69
## PACF 0.64 -0.11 -0.33 -0.38 -0.44 -0.29 -0.28 -0.07 -0.12 -0.03  0.03 -0.01
##      [,13] [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24]
## ACF   0.61  0.33  0.00 -0.37 -0.63 -0.72 -0.60 -0.33  0.05  0.35  0.63  0.72
## PACF  0.09 -0.05  0.01 -0.13 -0.19 -0.17 -0.15 -0.11 -0.03 -0.17 -0.04 -0.01
##      [,25] [,26] [,27] [,28] [,29] [,30] [,31] [,32] [,33] [,34] [,35] [,36]
## ACF   0.60  0.32 -0.02 -0.35 -0.59 -0.69 -0.58 -0.33  0.01  0.35  0.57  0.65
## PACF -0.02 -0.04 -0.04  0.02 -0.01 -0.03 -0.01 -0.06 -0.02 -0.01 -0.03 -0.03
##      [,37] [,38] [,39] [,40] [,41] [,42] [,43] [,44] [,45] [,46] [,47] [,48]
## ACF   0.57  0.30  0.00 -0.36 -0.58 -0.65 -0.55 -0.31  0.03  0.33  0.57  0.65
## PACF  0.03 -0.05  0.09 -0.08  0.00 -0.01 -0.01  0.00  0.00 -0.03  0.01  0.04
monthplot(ndvi.obs1, xlab="MESES", ylab="NDVI", type="l")

MODELOS

Podemos pensar em modelos autoregressivos de ordem até 6, ainda podemos considerar comportamento sazonal. Este não será considerado no momento. 

mod01 = arima(ndvi.obs1, order = c(6,0,0))
mod01
## 
## Call:
## arima(x = ndvi.obs1, order = c(6, 0, 0))
## 
## Coefficients:
##          ar1     ar2     ar3      ar4      ar5      ar6  intercept
##       0.2318  0.1027  0.0214  -0.0937  -0.3460  -0.3047     -1e-04
## s.e.  0.0532  0.0514  0.0517   0.0516   0.0514   0.0534      8e-04
## 
## sigma^2 estimated as 0.0003568:  log likelihood = 822.05,  aic = -1628.11
mod02 = arima(ndvi.obs1, order = c(6,0,1))
mod02
## 
## Call:
## arima(x = ndvi.obs1, order = c(6, 0, 1))
## 
## Coefficients:
##          ar1      ar2      ar3      ar4      ar5      ar6      ma1  intercept
##       0.7509  -0.0997  -0.0539  -0.0991  -0.2597  -0.0918  -0.6483     -1e-04
## s.e.  0.1027   0.0727   0.0694   0.0689   0.0686   0.0816   0.0909      4e-04
## 
## sigma^2 estimated as 0.0003247:  log likelihood = 836.96,  aic = -1655.93
mod03 = arima(ndvi.obs1, order = c(6,0,2))
mod03
## 
## Call:
## arima(x = ndvi.obs1, order = c(6, 0, 2))
## 
## Coefficients:
##          ar1     ar2      ar3      ar4      ar5      ar6      ma1      ma2
##       0.1393  0.5049  -0.1573  -0.1497  -0.3714  -0.1725  -0.0248  -0.5532
## s.e.  0.2783  0.2151   0.0732   0.0646   0.0622   0.1201   0.2760   0.1871
##       intercept
##          -1e-04
## s.e.      4e-04
## 
## sigma^2 estimated as 0.0003253:  log likelihood = 836.64,  aic = -1653.28

PREVISÃO

# horizonte de previsão
m = 6
preditos = predict(mod02, n.ahead = m)
par(mfrow = c(1,1), mar=c(4,3,1,1), mgp=c(1.6,.6,0), pch=19)
ts.plot(ndvi.obs1, preditos$pred, col=1:2, xlim = c(1995,2010), lwd=2, ylab="NDVI", xlab="", type="b")
U = preditos$pred+preditos$se; L = preditos$pred-preditos$se
xx = c(time(U), rev(time(U))); yy = c(L, rev(U))
polygon(xx, yy, border = 8, col = gray(.6, alpha = .2))
lines(preditos$pred, type="p", col=2)
abline(h=0)
grid()

Como proceder para fazer as previsões nos dados originais? em situações que nem esta não é possível aplicar a trnasformação inversa que utilizamos nos dados. Significa que este tipo de transformação somente interessa caso se queira estudar a série sem pretenções de previsão.