# Plot of Co2 data data(co2) plot(co2, ylab = expression("Atmospheric concentration of CO"[2]),las = 1) title(main = "co2 data set") yt <- as.vector(co2) nt <- length(yt) tt <-seq(1:nt) # Plots of lowess and first difference co2smooth<-yt-lowess(tt,yt,f=0.2)$y co2diff<-diff(yt) par(mfrow=c(2,2)) ts.plot(co2smooth) mtext("lowess, f=0.2") acf(co2smooth,lag=50) ts.plot(co2diff) mtext("First difference") acf(co2diff,lag=50) # # function to evaluate log-likelihood with approx lik<-function(w){x<-NULL ; ss<-sum(yt*yt); for(i in 1: k){wi<-w[i] ; r<-sum(yt*cos(wi*tt))^2+sum(yt*sin(wi*tt))^2; x<-c(x,log(1-2*r/(nt*ss))*(2-nt)/2) } return(x)} # # Plots of periodograms par(mfrow=c(3,2)) yt<-co2smooth lam<-seq(2,16,length=100) k<-length(w<-2*pi/lam) plam<-lik(w) plot(lam,plam,type='l',xlab='lambda',ylab='log lik') mtext("lowess") plot(lam,exp(plam-max(plam)),type='l',xlab='lambda',ylab='lik') mtext("lowess") # Removing cycle of periodicity 12 from the data. a<-sum(yt*cos(2*pi*tt/12))*2/nt ; b<-sum(yt*sin(2*pi*tt/12))*2/nt yt<-yt-(fit1<-a*cos(2*pi*tt/12) + b*sin(2*pi*tt/12)) plam<-lik(w) plot(lam,plam,type='l',xlab='lambda',ylab='log lik') mtext("Fitted cycle of period 12") plot(lam,exp(plam-max(plam)),type='l',xlab='lambda',ylab='lik') mtext("Fitted cycle of period 12") # # Removing a cycle of periodicity 6 a<-sum(yt*cos(2*pi*tt/6))*2/nt ; b<-sum(yt*sin(2*pi*tt/6))*2/nt yt<-yt-(fit2<-a*cos(2*pi*tt/6) + b*sin(2*pi*tt/6)) plam<-lik(w) plot(lam,plam,type='l',xlab='lambda',ylab='log lik') mtext("Fitted cycle of period 6") plot(lam,exp(plam-max(plam)),type='l',xlab='lambda',ylab='lik') mtext("Fitted cycle of period 6") # Data, estimated cycles and residuals par(mfrow=c(3,2)) yt<-co2smooth yl<-range(yt) ts.plot(yt,ylab='smoothed co2 data',ylim=yl) ts.plot(fit1,ylab='fitted cycle lambda=12',ylim=yl) ts.plot(fit2,ylab='fitted cycle lambda=6',ylim=yl) ts.plot(fit1+fit2,ylab='sum of fitted cycles',ylim=yl) ts.plot(residuals<-yt-(fit1+fit2),ylab="residuals",ylim=yl) acf(residuals,lag=50) # par(mfrow=c(3,2)) yt<-co2diff tt<-seq(1:(length(yt))) lam<-seq(2,16,length=100) k<-length(w<-2*pi/lam) plam<-lik(w) plot(lam,plam,type='l',xlab='lambda',ylab='log lik') mtext("diff") plot(lam,exp(plam-max(plam)),type='l',xlab='lambda',ylab='lik') mtext("diff") # a<-sum(yt*cos(2*pi*tt/12))*2/nt ; b<-sum(yt*sin(2*pi*tt/12))*2/nt yt<-yt-(fit1<-a*cos(2*pi*tt/12) + b*sin(2*pi*tt/12)) plam<-lik(w) plot(lam,plam,type='l',xlab='lambda',ylab='log lik') mtext("Fitted cycle of period 12") plot(lam,exp(plam-max(plam)),type='l',xlab='lambda',ylab='lik') mtext("Fitted cycle of period 12") # a<-sum(yt*cos(2*pi*tt/6))*2/nt ; b<-sum(yt*sin(2*pi*tt/6))*2/nt yt<-yt-(fit2<-a*cos(2*pi*tt/6) + b*sin(2*pi*tt/6)) plam<-lik(w) plot(lam,plam,type='l',xlab='lambda',ylab='log lik') mtext("Fitted cycle of period 6") plot(lam,exp(plam-max(plam)),type='l',xlab='lambda',ylab='lik') mtext("Fitted cycle of period 6") dev.off() # par(mfrow=c(3,2)) yt<-co2diff yl<-range(yt) ts.plot(yt,ylab='smoothed co2 data',ylim=yl) ts.plot(fit1,ylab='fitted cycle lambda=12',ylim=yl) ts.plot(fit2,ylab='fitted cycle lambda=6',ylim=yl) ts.plot(fit1+fit2,ylab='sum of fitted cycles',ylim=yl) ts.plot(residuals<-yt-(fit1+fit2),ylab="residuals",ylim=yl) acf(residuals,lag=50) # lik<-function(w){ x<-NULL; for (i in 1:k) { ss<-sum(yt*yt); wi<-w[i]; ci<-cos(wi*tt); si<-sin(wi*tt); C<-sum(ci*ci); S<-sum(si*si); E<-sum(ci*si); D<-C*S-E*E; cf<-sum(yt*ci); sf<-sum(yt*si); r<-(cf*cf*S-2*cf*sf*E+sf*sf*C)/D; x<-c(x,-log(D)/2+log(1-r/ss)*(2-nt)/2); }; x}; par(mfrow=c(3,2)) yt<-co2diff tt<-seq(1:(length(yt))) lam<-seq(2.5,16,length=100) k<-length(w<-2*pi/lam) plam<-lik(w) plot(lam,plam,type='l',xlab='lambda',ylab='log lik') mtext("diff") plot(lam,exp(plam-max(plam)),type='l',xlab='lambda',ylab='lik') mtext("diff") # a<-sum(yt*cos(2*pi*tt/12))*2/nt ; b<-sum(yt*sin(2*pi*tt/12))*2/nt yt<-yt-(fit1<-a*cos(2*pi*tt/12) + b*sin(2*pi*tt/12)) plam<-lik(w) plot(lam,plam,type='l',xlab='lambda',ylab='log lik') mtext("Fitted cycle of period 12") plot(lam,exp(plam-max(plam)),type='l',xlab='lambda',ylab='lik') mtext("Fitted cycle of period 12") # a<-sum(yt*cos(2*pi*tt/6))*2/nt ; b<-sum(yt*sin(2*pi*tt/6))*2/nt yt<-yt-(fit2<-a*cos(2*pi*tt/6) + b*sin(2*pi*tt/6)) plam<-lik(w) plot(lam,plam,type='l',xlab='lambda',ylab='log lik') mtext("Fitted cycle of period 6") plot(lam,exp(plam-max(plam)),type='l',xlab='lambda',ylab='lik') mtext("Fitted cycle of period 6") # par(mfrow=c(3,2)) yt<-co2diff yl<-range(yt) ts.plot(yt,ylab='smoothed co2 data',ylim=yl) ts.plot(fit1,ylab='fitted cycle lambda=12',ylim=yl) ts.plot(fit2,ylab='fitted cycle lambda=6',ylim=yl) ts.plot(fit1+fit2,ylab='sum of fitted cycles',ylim=yl) ts.plot(residuals<-yt-(fit1+fit2),ylab="residuals",ylim=yl) acf(residuals,lag=50) # par(mfrow=c(4,2)) per<-spec.pgram(co2diff,taper=0,pad=0,detrend=F,demean=F) lam<-1/per$freq plam<-per$spec i<-2