# co2 analysis
# no need to tell how much to skip 
# "comment lines are ignored in the file '#'
#rm(list=ls())
buffer<-read.table("ftp://aftp.cmdl.noaa.gov/products/trends/co2/co2_mm_mlo.txt")
time<-buffer[,3]
co2<-buffer[,5]
n<-length(time)
co2m<-mean(co2)
par(mfrow=c(1,1))
plot(time,co2,typ='l',xlab="time",ylab='CO2 [ppm]')


############################################
# harmonic analysis
############################################
# fitting the amplitude of the annual cycle
# with linear regression
# NOTE: sine and cosine functions are
# orthogonal functions i.e. their time series 
# 'appear' uncorrelated 
# making it possible to use multiple linear regression
f1<-1/12
f2<-1/6
f3<-1/4
index<-1:n
s1<-sin(2*pi*f1*index)
c1<-cos(2*pi*f1*index)
#s2<-sin(2*pi*f2*index)
#c2<-cos(2*pi*f2*index)
#s3<-sin(2*pi*f3*index)
#c3<-cos(2*pi*f3*index)
#time<-buffer[,3]

res<-lm(buffer[,5] ~ s1 + c1  )


print(summary(res))
lines(buffer[,3],res$fitted,col=2)

lines(buffer[,3],co2m+res$resid,col=3)

# detrending data via first-order differences
# calculated the difference between neighboring values
# 
d1<-res$resid[2:n]-res$resid[1:(n-1)]
lines(time[2:n],d1[2:n]+co2m,col='blue')

# use the acf() function to study the autocorrelation 
# in the residuals


# use the running mean function to remove the annual cycle
# then the trend