##设置工作目录
##setwd()
##设置图片保存地址
#fig.dir = "" 
data(AirPassengers)
AP <- AirPassengers
pdf(paste(fig.dir,"AirPassengers.pdf",sep=""))
plot(AP, ylab = "Passengers (1000's)")
dev.off()




####Maine Unemployment data
www <- "http://www.math.pku.edu.cn/teachers/xirb/Courses/TimeSeries/IntroTimeSeriesRData/Maine.dat"
Maine.month <- read.table(www, header = TRUE)
attach(Maine.month)
Maine.month.ts <- ts(unemploy, start = c(1996, 1), freq = 12)
Maine.annual.ts <- aggregate(Maine.month.ts)/12
pdf(paste(fig.dir,"MainUnemployment.pdf",sep=""))
par(mfrow=c(2,1))
plot(Maine.month.ts, ylab = "unemployed (%)",main="Monthly")
plot(Maine.annual.ts, ylab = "unemployed (%)",main="Annually")
dev.off()
detach(Maine.month)


z= decompose(Maine.month.ts)
plot(z$trend)
plot(z$seasonal,las=2)



#### global warming data
www <- "http://www.math.pku.edu.cn/teachers/xirb/Courses/TimeSeries/IntroTimeSeriesRData/global.dat"
Global <- scan(www)
Global.ts <- ts(Global, st = c(1856, 1), end = c(2005, 12),
fr = 12)
Global.annual <- aggregate(Global.ts, FUN = mean)
pdf(paste(fig.dir,"GobalTemperature.pdf",sep=""))
par(mfrow=c(2,1))
plot(Global.ts,ylab="Monthly Mean Temperature")
plot(Global.annual,ylab="Annually Mean Temperature")
dev.off()

set.seed(1346)
x = arima.sim(n = 1000, list(ar = c(0.98),sd = sqrt(0.01)))
plot(x,ylab="Value")
ind = c(120:260)
pdf(paste(fig.dir,"GobalTemperaturSimData_partial.pdf",sep=""))
plot(ts(x[ind]),ylab="Value")
dev.off()
pdf(paste(fig.dir,"GobalTemperaturSimData_all.pdf",sep=""))
plot(x,ylab="Value")
lines(ts(x[ind]),x=c(ind),col="red")
abline(v=ind[1])
abline(v=ind[length(ind)])
dev.off()



##############Multiple time series###############
www <- "http://www.math.pku.edu.cn/teachers/xirb/Courses/TimeSeries/IntroTimeSeriesRData/cbe.dat"
CBE <- read.table(www, header = T)
Elec.ts <- ts(CBE[, 3], start = 1958, freq = 12)
Beer.ts <- ts(CBE[, 2], start = 1958, freq = 12)
Choc.ts <- ts(CBE[, 1], start = 1958, freq = 12)
pdf(paste(fig.dir,"AustrilianChoBeerElec.pdf",sep=""))
plot(cbind(Elec.ts, Beer.ts, Choc.ts),main="Chocolate, Beer, and Electricity Production: 1958-1990")
dev.off()

pdf(paste(fig.dir,"AustrilianElec_AP.pdf",sep=""))
AP.elec <- ts.intersect(AP, Elec.ts)
AP <- AP.elec[,1]; Elec <- AP.elec[,2]
layout(1:3)
plot(AP, main = "", ylab = "Air passengers / 1000's")
plot(Elec, main = "", ylab = "Electricity production / MkWh")
plot(as.vector(AP), as.vector(Elec),
xlab = "Air passengers / 1000's",
ylab = "Electricity production / MWh")
abline(reg = lm(Elec ~ AP))
dev.off()



####### now the data from "应用时间序列分析" ############
source('http://www.math.pku.edu.cn/teachers/xirb/Courses/TimeSeries/code/data.r')
source("http://www.math.pku.edu.cn/teachers/xirb/Courses/TimeSeries/code/atsa01.r")

demo.coal()


######## #############
data(AirPassengers)
AP <- AirPassengers
AP = log(AP)
### Use decompose
par(mfrow=c(2,1),lwd=2)
plot(AP)
z = decompose(AP)
lines(z$trend,col="blue")
lines(z$trend+z$seasonal,col="red")
legend("topleft",lty=1,col=c("black","blue","red"),legend=c("Observed Data","Trend","Fitted"))
plot(z$seasonal,col="blue",ylim=c(-0.03,0.03))
lines(z$random,col="green")
legend("topleft",lty=1,col=c("black","blue"),legend=c("Random","Seasonal"))

### Use stl
z = stl(AP,s.window="periodic",s.degree=0)
par(mfrow=c(2,1),lwd=2)
plot(AP)
lines(z$time.series[,2],col="blue")
lines(z$time.series[,2]+z$time.series[,1],col="red")
legend("topleft",lty=1,col=c("black","blue","red"),legend=c("Observed Data","Trend","Fitted"))
plot(z$time.series[,1],col="red")
lines(z$time.series[,3],col="green")
legend("topleft",lty=1,col=c("black","blue"),legend=c("Random","Seasonal"))