soi=read.table("datos/soi.dat", dec=".",header=F)$V1 rec=read.table("datos/recruit.dat",dec=".",header=F)$V1 soi.ts<- ts(soi, start=c(1950, 1), frequency=12) rec.ts <- ts(rec, start=c(1950, 1), frequency=12) par(mfrow=c(2,1),mar=c(3, 4, 3, 2)) plot(soi.ts,main="Southern Oscillation Index (El Niño)",ylab="") plot(rec.ts,main="Pacific Recruitment",ylab="") ## Ruido Blanco, con la misma longitud y sd que los datos "soi" wn=rnorm(length(soi),0,sd(soi)) wn.ts=ts(wn,start=c(1950,1), frequency=12) plot(soi.ts,main="Southern Oscillation Index (El Niño)",ylab="") plot(wn.ts,main="White Noise",ylab="") ## Proceso MA (Moving Average) mawn=stats::filter(wn.ts,sides=2,rep(1,3)/3) plot(wn.ts,main="White Noise",ylab="") plot(mawn,main="Moving Average Process",ylab="") ## Proceso AR (AutoRegresivo) w0=c(rnorm(50,0,sd(soi)),wn) x=stats::filter(w0,filter=c(1,-0.9),method="recursive") x.ts=ts(x[-(1:50)],start=c(1950,1), frequency=12) par(mfrow=c(3,1),mar=c(3, 4, 3, 2)) plot(wn.ts,main="Ruido Blanco",ylab="") plot(mawn,main="Proceso MA (Media Móvil)",ylab="") plot(x.ts,main="Proceso AR (Autoregresivo)") ## Proceso no estacionario (Señal + ruido) t=1:240 y=2*cos(2*pi*t/50+0.6*pi) wn=rnorm(240,0,1) par(mfrow=c(3,1)) plot.ts(y) plot.ts(y+wn) plot.ts(y+5*wn) x = ts(cbind(soi,rec)) s = spec.pgram(x, kernel("daniell",9), taper=0) s$df # df = 35.8625 f = qf(.999, 2, s$df-2) # f = 8.529792 c = f/(18+f) # c = 0.3188779 plot(s, plot.type = "coh", ci.lty = 2) abline(h = c) ## ACF y PACF wn=rnorm(240,0,1) par(mfrow=c(2,1)) acf(wn,main="ACF ruido blanco") pacf(wn,main="PACF ruido blanco") acf(mawn,main="ACF Proceso MA",na.action=na.exclude) pacf(mawn,main="PACF Proceso MA",na.action=na.exclude) acf(x.ts,main="ACF Proceso AR") pacf(x.ts,main="PACF Proceso AR") acf(y+wn,main="ACF Proceso Señal+ruido") pacf(y+wn,main="PACF Proceso Señal+ruido") acf(y+5*wn,main="ACF Proceso Señal+ruido") pacf(y+5*wn,main="PACF Proceso Señal+ruido") ## ccf entre ruidos blancos t=1:240 w1=rnorm(240,0,1) w2=rnorm(240,0,1) ccf(w1,w2,ylab="ccf",main="Correlación cruzada") ## ccf entre señales de igual frecuencia t=1:240 y=2*cos(2*pi*t/50+0.6*pi)+rnorm(240,0,1) z=2*cos(2*pi*t/50+0.6*pi)+rnorm(240,0,1) ccf(y,z,ylab="ccf",main="Correlación cruzada") ## ccf entre soi y reclutamiento par(mfrow=c(3,1)) acf(soi, 50) acf(rec, 50) ccf(soi, rec, 50) require(TSA) par(mfrow=c(3,1)) acf(soi, 30) pacf(soi, 30) t=1:length(soi) soi.det=residuals(lm(soi~t)) acf(soi.det, 30) pacf(soi.det, 30) pw=prewhiten(soi.det,rec,x.model=arima(soi,order=c(1,0,0)), ylab="CCF") pw=prewhiten(soi,rec,x.model=arima(soi,order=c(1,0,0)), ylab="CCF") prewhiten(rec,soi.det,ylab="CCF") acf(rec.ts) pacf(rec.ts) mrec=arima(rec.ts,order=c(2,0,0)) tsdiag(mrec) rmr=residuals(mrec) pw=prewhiten(soi.ts,rmr,x.model=arima(soi,order=c(1,0,0)), ylab="CCF") ## Regresion (m1=arima(soi.det,order=c(1,0,0))) tsdiag(m1) lr=length(rec) dfm2=lm(rec[-(1:5)]~rec[-c(1:4,lr)]+soi[-c((lr-4):lr)]) summary(dfm2) require(dynlm) # Idem con dynlm. Hay que declarar las variables como ts, si no no funciona dfm <- dynlm(rec.ts ~ L(rec.ts, 1)+L(rec.ts, 2)+L(soi.ts,5)) summary(dfm) plot(residuals(dfm)) acf(residuals(dfm)) pacf(residuals(dfm)) mr=arima(residuals(dfm2),order=c(2,0,0)) tsdiag(mr) # Mixed Effects Models and Extensions in Ecology with R (2009) # Zuur, Ieno, Walker, Saveliev and Smith. Springer # This file was produced by Alain Zuur (highstat@highstat.com) # www.highstat.com # This program is free software; you can redistribute it and/or modify # it under the terms of the GNU General Public License as published by # the Free Software Foundation; either version 2 of the License, or # (at your option) any later version. # # This program is distributed in the hope that it will be useful, # but WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # GNU General Public License for more details. ## Lectura de los datos y representación gráfica Hawaii<-read.table("http://dl.dropbox.com/u/7610774/Datos/zuur/Hawaii.txt",header=T) Hawaii$Birds<-sqrt(Hawaii$Moorhen.Kauai) plot(Hawaii$Year,Hawaii$Birds,xlab="Year", ylab="Moorhen abundance on Kauai",type="b") ## Ajuste modelo lineal library(nlme) M0 <- gls(Birds ~ Rainsoi=read.table("datos/soi.dat", dec=".",header=F)$V1 rec=read.table("datos/recruit.dat",dec=".",header=F)$V1 soi.ts<- ts(soi, start=c(1950, 1), frequency=12) rec.ts <- ts(rec, start=c(1950, 1), frequency=12) par(mfrow=c(2,1),mar=c(3, 4, 3, 2)) plot(soi.ts,main="Southern Oscillation Index (El Niño)",ylab="") plot(rec.ts,main="Pacific Recruitment",ylab="") ## Ruido Blanco, con la misma longitud y sd que los datos "soi" wn=rnorm(length(soi),0,sd(soi)) wn.ts=ts(wn,start=c(1950,1), frequency=12) plot(soi.ts,main="Southern Oscillation Index (El Niño)",ylab="") plot(wn.ts,main="White Noise",ylab="") ## Proceso MA (Moving Average) mawn=filter(wn.ts,sides=2,rep(1,3)/3) plot(wn.ts,main="White Noise",ylab="") plot(mawn,main="Moving Average Process",ylab="") ## Proceso AR (AutoRegresivo) w0=c(rnorm(50,0,sd(soi)),wn) x=filter(w0,filter=c(1,-0.9),method="recursive") x.ts=ts(x[-(1:50)],start=c(1950,1), frequency=12) par(mfrow=c(3,1),mar=c(3, 4, 3, 2)) plot(wn.ts,main="Ruido Blanco",ylab="") plot(mawn,main="Proceso MA (Media Móvil)",ylab="") plot(x.ts,main="Proceso AR (Autoregresivo)") ## Proceso no estacionario (Señal + ruido) t=1:240 y=2*cos(2*pi*t/50+0.6*pi) wn=rnorm(240,0,1) par(mfrow=c(3,1)) plot.ts(y) plot.ts(y+wn) plot.ts(y+5*wn) x = ts(cbind(soi,rec)) s = spec.pgram(x, kernel("daniell",9), taper=0) s$df # df = 35.8625 f = qf(.999, 2, s$df-2) # f = 8.529792 c = f/(18+f) # c = 0.3188779 plot(s, plot.type = "coh", ci.lty = 2) abline(h = c) ## ACF y PACF wn=rnorm(240,0,1) par(mfrow=c(2,1)) acf(wn,main="ACF ruido blanco") pacf(wn,main="PACF ruido blanco") acf(mawn,main="ACF Proceso MA",na.action=na.exclude) pacf(mawn,main="PACF Proceso MA",na.action=na.exclude) acf(x.ts,main="ACF Proceso AR") pacf(x.ts,main="PACF Proceso AR") acf(y+wn,main="ACF Proceso Señal+ruido") pacf(y+wn,main="PACF Proceso Señal+ruido") acf(y+5*wn,main="ACF Proceso Señal+ruido") pacf(y+5*wn,main="PACF Proceso Señal+ruido") ## ccf entre ruidos blancos t=1:240 w1=rnorm(240,0,1) w2=rnorm(240,0,1) ccf(w1,w2,ylab="ccf",main="Correlación cruzada") ## ccf entre señales de igual frecuencia t=1:240 fall + Year, na.action = na.omit, data = Hawaii) summary(M0) ## ACF de los residuos E<-residuals(M0,type="normalized") ## OJO! los siguientes comandos permiten reintroducir los valores perdidos en la serie de los ## residuos; de no hacerse así, esta serie sería contínua y no correspondería a la secuencia ## real de los residuos. I<-!is.na(Hawaii$Birds) Efull<-vector(length=length(Hawaii$Birds)) Efull<-NA Efull[I]<-E acf(Efull,na.action = na.pass, main="Auto-correlation plot for residuals") ## Los p-valores del modelo no son fiables porque los residuos no son independientes. ## El siguiente modelo estima considerando que los residuos siguen un modelo AR1 M2<-gls(Birds ~ Rainfall + Year, na.action = na.omit, correlation = corAR1(form =~ Year), data = Hawaii) summary(M2) ## De modo equivalente: birds=ts(Hawaii$Birds, start=1956) rain=ts(Hawaii$Rainfall,start=1956) year=Hawaii$Year Mts1=arima(birds,xreg=data.frame(year,rain),order=c(1,0,0)) Mts1 cs1 <- corARMA(c(0.2), p = 1, q = 0) cs2 <- corARMA(c(0.3, -0.3), p = 2, q = 0) M3arma1<-gls(Birds~Rainfall+Year,na.action=na.omit, correlation=cs1,data = Hawaii) M3arma2<-gls(Birds~Rainfall+Year,na.action=na.omit, correlation=cs2,data = Hawaii) AIC(M3arma1,M3arma2)