I am working on creating a distributed lag model between monthly NDVI and rainfall values. The data contains a stack of monthly values from 1982 - 2013, and I have to create the lag model for any one pixel in the dataset (in this case [170,200,] (see code block). I have been trying to remove the trend and seasonality from the data but it hasn't worked thus far.
filename<-"E:/Semester 2/TSA/7/ndvi_1982_2013_sub"
filename2<-"E:/Semester 2/TSA/7/rain_1982_2013_sub"
ndvi<- read.ENVI(filename, headerfile=paste(filename, ".hdr", sep=""))
rain<- read.ENVI(filename2, headerfile=paste(filename, ".hdr", sep=""))
datan <- ts(ndvi[170,200,],start=c(1982,1),frequency=12)
datar <- ts(rain[170,200,],start=c(1982,1),frequency=12)
model <- lm(datan ~ datar)
##here i model the original data just to look at
##the ACF and PACF structure
par(mfrow=c(1,2))
acf(model$residuals,main="")
acf(model$residuals,type="p",main="")
ACF and PACF of the first regression model, clear seasonality, so I used stl to remove the seasonal and trend components
decompdatan <- stl(datan,"periodic")
decompdatar <- stl(datar,"periodic")
acf(decompdatan[[1]][,3],main="")
acf(decompdatan[[1]][,3],type='p',main="")
Here the ACF and PACF of the residuals. Although the seasonal component has been removed there is still an apparent seasonal flucatuation.
According to the assignment, I am supposed to move on to the pre-whitening stage after this step in order to determine the number of lags I should incorporate into the model, but the data is clearly not ready for this stage yet since there are some many significantly correlated lags still apparent in the data. Also, when I tried to iterate the number of lags within the arima model (for the pre-whitening), the AIC kept improving well past the orders of 8 and 9, but I have been told that the AIC should normally decrease until it reaches a sort of plateau, after which AIC values should again increase, thus making a clear choice of which order is most suitable.
I have tried other ways of removing the seasonal component such as using a different function and trying to create my own seasonal component using Fourier polynomials, but nothing seems to work. Is there something glaring I am doing wrong here?
Thanks.