rm(list=objects()) library(magrittr) ####################################################################################################################### ##################################exercice 1 ####################################################################################################################### sigma<-1/2 phi<-0.9 n<-1000 eps<-rnorm(n,0,sd=sigma) y<-arima.sim(n = n, list(ar = c(phi)),innov=eps) var(y) sigma^2/(1-phi^2) ?arima.sim plot(y,type='l') Nsimu <- 100 v <- array(0, dim=Nsimu) for(i in c(1:Nsimu)) { eps<-rnorm(n,0,sd=sigma) y<-arima.sim(n = n, list(ar = c(phi)),innov=eps) v[i] <- var(y) } boxplot(v) abline(h=sigma^2/(1-phi^2), col='red') ##################autocorrelation empirique et theorique lmax <- 40 rho_est<-acf(y,lag.max=lmax, plot = FALSE) rho_theo<-phi^c(0:lmax) rho_est rho_theo plot(rho_est$acf,type='b',pch=20) lines(rho_theo,col='red') ##################autocovariance empirique et th?orique gamma_est<-rho_est$acf*var(y) gamma_theo<-sigma^2/(1-phi^2)*phi^c(0:50) phi=0.8 sigma=1 gamma_theo<-sigma^2/(1-phi^2)*phi^c(0:50) gamma_theo gamma_est-gamma_theo plot(gamma_est,type='b',pch=20) lines(gamma_theo,col='red') ######estimateur de phi: phi_hat<-rho_est$acf[2] phi_hat ######intervalle de confiance sigma<-1/2 phi<-0.7 n<-500 estim_phi<-function(n,sigma,phi) { eps<-rnorm(n,0,sd=sigma) y<-arima.sim(n = n, list(ar = c(phi)),innov=eps) rho_est<-acf(y,lag.max=1, plot=FALSE) phi_hat<-rho_est$acf[2] return(phi_hat) } Nsimu<-500 phi_hat<-lapply(rep(n,Nsimu),estim_phi,sigma=sigma,phi=phi) phi_hat<-unlist(phi_hat) ###histogramme des estimations de phi hist(phi_hat,breaks=30) abline(v=phi,col='red') ###intervalle de confiances (2 approches) alpha<-0.05 #solution1 a<-sort(phi_hat)[floor(alpha/2*Nsimu)] b<-sort(phi_hat)[floor((1-alpha/2)*Nsimu)] c(a,b) #solution2 quantile(phi_hat,c(alpha/2,1-alpha/2)) Nsimu<-100 s_n<-seq(50,500,by=50) IC<-NULL for(n in s_n) { phi_h<-lapply(rep(n,Nsimu), estim_phi,sigma=sigma,phi=phi) phi_h<-unlist(phi_h) IC<-rbind(IC,quantile(phi_h,c(alpha/2, 1-alpha/2))) print(n) } plot(s_n,IC[,1],pch='-',ylim=range(IC),cex=3) points(s_n,IC[,2],pch='-',cex=3) l<-abs(IC[,2]-IC[,1]) plot(s_n,l,type='b',pch=20) conv<-1/sqrt(s_n) reg<-lm(l~conv-1) lines(s_n,reg$coeff/sqrt(s_n),col='red') ####################################################################################################################### ##################################exercice 2 ####################################################################################################################### # ####donn?es simul?es ainsi # sigma<-1/2 # n<-100 # eps<-rnorm(n,0,sd=sigma) # t<-c(1:n) # s<-cos(2*pi*t/10)+2*cos(4*pi*t/10)+2*cos(8*pi*t/10) # y<-s+eps # #y<-s+eps # plot(y,type='l') # spectrum(y) # abline(v=c(1,2,4)/10,col='red') # # # # # exercice2<-data.frame(t,y) # #setwd("/Users/yannig/Documents/Enseignement/2017-2018/Serie_temp/TP/TP5_data/exercice2.txt") # write.table(exercice2,file='exercice2.txt',col.names=T,row.names=F,quote=F,sep=';') ######################################################## #exercice2.txt #setwd("/Users/yannig/Documents/Enseignement/2017-2018/Serie_temp/TP/TP5_data/") data<-read.table("TP5_data/exercice2.txt",sep=';', header=T) par(mfrow=c(1,1)) plot(data$y,type='l') s<-spectrum(data$y) plot(s$freq,s$spec,type='l') ############construction de la serie de fourier a1<-which.max(s$spec) f1<-s$freq[a1] 1/f1 a2<-which.max(s$spec[-a1]) f2<-s$freq[a2] a3<-which.max(s$spec[-c(a1,a2)]) f3<-s$freq[a3] 1/c(f1,f2,f3) t <- c(1:length(data$y)) X1cos<-cos(2*pi*f1*data$t) X1sin<-sin(2*pi*f1*data$t) X2cos<-cos(2*pi*f2*data$t) X2sin<-sin(2*pi*f2*data$t) X3cos<-cos(2*pi*f3*data$t) X3sin<-sin(2*pi*f3*data$t) model1<-lm(data$y~X1cos+X1sin+X2cos+X2sin+X3cos+X3sin-1) summary(model1) plot(data$y,type='l') lines(model1$fitted, col='red') data0 <- data.frame(y=data$y, X1cos=X1cos, X1sin=X1sin, X2cos=X2cos, X2sin=X2sin, X3cos=X3cos, X3sin=X3sin) model1<-lm(y~X1cos+X1sin+X2cos+X2sin+X3cos+X3sin-1, data=data0) summary(model1) predict(model1, newdata=data1) model2<-lm(y~X1cos+X2cos+X3cos-1, data=data0) summary(model2) plot(data$y,type='l') lines(model2$fitted, col='red') t <- c((length(data$y)+1) : (length(data$y)+50) ) X1cos<-cos(2*pi*f1*t) X2cos<-cos(2*pi*f2*t) X3cos<-cos(2*pi*f3*t) data1 <- data.frame(X1cos=X1cos, X2cos=X2cos, X3cos=X3cos) pred <- predict(model2, newdata=data1) length(pred) plot(t, pred, type='l') ####################################################################################################################### ##################################exercice 3 ####################################################################################################################### set.seed(110) n<-500 sigma<-1 eps<-rnorm(n,0,sd=sigma) x1<-arima.sim(n = n, list(ar = c(5/6,-3/4)),innov=eps) plot(x1,type='l') par(mfrow=c(1,2)) acf(x1) pacf(x1) n<-500 sigma<-1 eps<-rnorm(n,0,sd=sigma) x2<-arima.sim(n = n, list(ar = c(0.7),innov=eps)) plot(x2,type='l') par(mfrow=c(1,2)) acf(x2) pacf(x2) n<-500 sigma<-1 eps<-rnorm(n,0,sd=sigma) x3<-arima.sim(n = n, list(ma = c(0.8,0.5,-3,1,4,5),innov=eps)) plot(x3,type='l') par(mfrow=c(1,2)) acf(x3) pacf(x3) n<-500 sigma<-2 eps<-rnorm(n,0,sd=sigma) x4<-arima.sim(n = n, list(ma = c(1:20)/100,innov=eps)) plot(x4,type='l') par(mfrow=c(1,2)) acf(x4,lag.max=30) pacf(x4,lag.max=30) exercice3<-cbind(x1,x2,x3,x4) write.table(exercice3,file='exercice3.txt',col.names=T,row.names=F,quote=F,sep=';') ########################################################################################## #####corrige ########################################################################################## setwd("/Users/yannig/Documents/Enseignement/2017-2018/Serie_temp/TP/TP5_data/") data<-read.table('exercice3.txt',header=T,sep=';') data[1,] attach(data) ###########x1 x1<-ts(x1) plot(x1) par(mfrow=c(1,2)) acf(x1) pacf(x1) mean(x1) sd(x1) lm(x1~1)%>%summary() ?arima x1.model<-arima(x1, order = c(2,0,0), method = c("ML"), SSinit = c("Rossignol2011"), optim.method = "BFGS", include.mean = F) x1.model_3<-arima(x1, order = c(3,0,0), method = c("ML"), SSinit = c("Rossignol2011"), optim.method = "BFGS", include.mean = F) x1.model_5<-arima(x1, order = c(5,0,0), method = c("ML"), SSinit = c("Rossignol2011"), optim.method = "BFGS", include.mean = F) x1.model$aic x1.model_3$aic x1.model_5$aic #log likelihood = -695.04, aic = 1396.07 #log likelihood = -691.62, aic = 1395.24 x1.model summary(x1.model) horizon<-10 par(mfrow=c(1,1)) x1.forecast<-predict(x1.model,n.ahead = horizon,se.fit =F) plot(x1,xlim=c(1,nrow(data)+horizon)) lines(nrow(data)*c(1:horizon),x1.forecast,col='red') names(x1.model) x1.model$aic -2*x1.model$loglik+2*(2+1) # phi<-ARMAtoMA(ar = x1.model$coef, ma=0, 12) # plot(phi,type='l') # rho<-acf(x1,lag.max=2)$acf # rho # gamma_0<-x1.model$sigma2/(1-rho[-1]%*%x1.model$coef) # gamma_0 ###########x2 x2<-ts(x2) par(mfrow=c(1,1)) plot(x2) spectrum(x2) par(mfrow=c(1,2)) acf(x2) pacf(x2) par(mfrow=c(1,1)) x2.model<-arima(x2, order = c(1,0,0), method = c("ML"),SSinit = c("Rossignol2011"),optim.method = "BFGS",include.mean = F) x2.model horizon<-10 x2.forecast<-predict(x2.model,n.ahead = horizon,se.fit =F) par(mfrow=c(1,1)) plot(x2,xlim=c(1,nrow(data)+horizon)) lines(nrow(data)*c(1:horizon),x2.forecast,col='red') phi<-ARMAtoMA(ar = x2.model$coef, ma=0, 12) phi 0.7^c(1:12) plot(phi,type='l') ###########x3 x3<-ts(x3) par(mfrow=c(1,2)) acf(x3) pacf(x3) mean(x3) sd(x3) t.test(x3,alternative=c("two.sided")) x3.model<-arima(x3, order = c(0,0,5), method = c("ML"), SSinit = c("Rossignol2011"), optim.method = "BFGS",include.mean = F) #log likelihood = -1535.67, aic = 3083.35 #log likelihood = -1537.47, aic = 3082.95 x3.model horizon<-10 par(mfrow=c(1,1)) x3.forecast<-predict(x3.model,n.ahead = horizon,se.fit =F) plot(x3,xlim=c(1,nrow(data)+horizon)) lines(nrow(data)*c(1:horizon),x3.forecast,col='red') ###########x4 x4<-ts(x4) acf(x4, lag.max=30) pacf(x4, lag.max=30) mean(x4) x4.model<-arima(x4, order = c(0,0,20), method = c("ML"),SSinit = c("Rossignol2011"),optim.method = "BFGS",include.mean = F) x4.model horizon<-10 x4.forecast<-predict(x4.model,n.ahead = horizon,se.fit =F) plot(x4,xlim=c(1,nrow(data)+horizon)) lines(nrow(data)*c(1:horizon),x4.forecast,col='red') # #################################################################################################### # ##################################################cas non stationnaire # #################################################################################################### # # estim_phi2<-function(n,sigma) # { # eps<-rnorm(n,0,sd=sigma) # a<-1 # puiss<-function(a,k){c(a^c(k:1),rep(0,n-k))} # M<-lapply(c(1:n),puiss,a=a)%>%unlist%>%matrix(nrow=n,ncol=n,byrow=T) # y<-M%*%eps # rho_est<-acf(y,lag.max=1) # phi_hat<-rho_est$acf[2] # return(phi_hat) # } # # sigma<-1/2 # n<-500 # alpha<-0.05 # Nsimu<-500 # phi_hat<-lapply(rep(n,Nsimu),estim_phi2,sigma=sigma) # phi_hat<-unlist(phi_hat) # # ###histogramme des estimations de phi # hist(phi_hat,breaks=30) # abline(v=phi,col='red') # # # quantile(phi_hat,c(alpha/2,1-alpha/2)) # # # # # # # # # # # # # s_n<-seq(50,500,by=50) # IC<-NULL # for(n in s_n) # { # phi_h<-lapply(rep(n,Nsimu),estim_phi2,sigma=sigma) # phi_h<-unlist(phi_h) # IC<-rbind(IC,quantile(phi_h,c(alpha/2,1-alpha/2))) # print(n) # } # # # plot(s_n,IC[,1],pch='-',ylim=range(IC),cex=3) # points(s_n,IC[,2],pch='-',cex=3) # # # l<-abs(IC[,2]-IC[,1]) # plot(s_n,l,type='b',pch=20) # conv<-1/sqrt(s_n) # reg<-lm(l~conv-1) # lines(s_n,reg$coeff/sqrt(s_n),col='red') spectrum(y,kernel("daniell", c(3, 3))) par(mfrow=c(2,2)) plot(y07,type='l',ylim=range(y03,y09)) spectrum(y03,kernel("daniell", c(3, 3))) plot(y09,type='l') spectrum(y09,kernel("daniell", c(3, 3)))