第一题

1.安装时间序列包TSA, zoo, xts, quantmod, tseries, lubridate, dplyr, forecast, TTR

install.packages("TSA")
install.packages("zoo")
install.packages("xts")
install.packages("quantmod")
install.packages("tseries")
install.packages("lubridate")
install.packages("dplyr")
install.packages("forecast")
install.packages("TTR")
library("TSA")
## 
## 载入程辑包：'TSA'
## The following objects are masked from 'package:stats':
## 
##     acf, arima
## The following object is masked from 'package:utils':
## 
##     tar
library("zoo")
## 
## 载入程辑包：'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
library("xts")
library("quantmod")
## 载入需要的程辑包：TTR
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
library("tseries")
library("lubridate")
## 
## 载入程辑包：'lubridate'
## The following objects are masked from 'package:base':
## 
##     date, intersect, setdiff, union
library("dplyr")
## 
## 载入程辑包：'dplyr'
## The following objects are masked from 'package:xts':
## 
##     first, last
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library("forecast")
## Registered S3 methods overwritten by 'forecast':
##   method       from
##   fitted.Arima TSA 
##   plot.Arima   TSA
library("TTR")

第二题

  1. 模拟趋势平稳过程\(x_t=\mu+\beta t+w_t,\quad w_t\sim N(0,1)\)

  1. 取\(\beta=-0.3\)，在\(\mu=0\)和\(\mu=3\)时分别画出时序图（同一张图上）相比较

  2. 取\(\beta=-0.4\)，\(\mu=0.5\)，绘制时序图及自相关函数图；对该序列进行一阶差分，绘制差分后序列的时序图及自相关函数图，比较两种数据下的图像

#(1)
set_x <-function(mu,beta,n,seed){
  set.seed(seed)
  w <- rnorm(n,0,1)
  x <- numeric(n)
  for (t in 1:length(w)) x[t] <- mu+beta*t+w[t]
  return(x)
}

x1 <- set_x(mu=0,beta=-0.3,n=100,seed=1234)
x2 <- set_x(mu=3,beta=-0.3,n=100,seed=1234)

plot(x1, type="l", xlab="time", ylab="",col="red",main=bquote(X[t] ~ '='~ mu ~'+' ~ beta ~'*'~t~'+'~w[t] ~ ' with n = 100'),ylim=c(-40,5))
lines(x2,type="l", col="blue")
legend("bottomleft",col=c("red", "blue"), lty=c(1,1), lwd=c(2,2), legend=c("mu=0", "mu=3"))

#(2)
x <- set_x(mu=0.5,beta=-0.4,n=100,seed=1234)
par(mfrow=c(2,2))
plot(x, type="l", xlab="time", ylab="",col="red",main=bquote(X[t] ~ '='~ mu ~'+' ~ beta ~'*'~t~'+'~w[t] ~ 'mu=0.5'),ylim=c(-40,5))
acf(x)

x_diff <- diff(x)
plot(x_diff, type="l", xlab="time", ylab="",col="red",main="diff(x)")
acf(x_diff)

第三题

  1. 北京市1995年-2000年的月平均气温数据如附件data1.csv所示
  1. 绘制时序图，直观考察序列的确定性特征
#(1)
data1<-read.csv("data1.csv",header=TRUE,sep=",")
colnames(data1) <- c("month","1995","1996","1997","1998","1999","2000")

data1_ts<-as.matrix(data1[,-1])
data1_ts<-as.vector(data1_ts)
x.tem<-ts(data1_ts,start=c(1995,1),frequency = 12)
plot(x.tem,ylab=paste("Monthly temperature"),type="o",pch=16,nc=1,main="北京市1995年-2000年的月平均气温")

monthplot(x.tem, xlab="", ylab="")

seasonplot(x.tem, year.labels="TRUE", main="")

从时序图来看，数据具有明显的季节性，但趋势性不明显。不同年的同一个月平均气温相差不大，有不明显的下降趋势，其中7月气温最高，1月气温最低，从月度图和季节图也可以看出这点。

  1. 比照时间序列分解的例题采用不同的方法拟合趋势和季节性

滑动平均法拟合（decompose函数）

#(2)
#平滑处理
opar <- par(no.readonly = TRUE)
par(mfrow=c(2,2))
ylim <- c(min(x.tem)-1,max(x.tem)+1)
# 趋势图
plot(x.tem,main="Raw time series")
#ma
plot(ma(x.tem,10),main="Simple Moving Averages (k=10)",ylab="monthly temperature",ylim=ylim)
#SMA
plot(SMA(x.tem,n=15),ylab="monthly temperature",main="Simple Moving Averages (k=15)")
#rollmean
plot(rollmean(x.tem,20),ylab="monthly temperature",main="Simple Moving Averages (k=20)",ylim=ylim)

par(opar)

#decompose和stl函数拟合趋势性和季节性(滑动平均去趋势)
par(mfrow=c(1,2))
d1=decompose(x.tem)
plot(d1)

d2=stl(x.tem,'periodic')
plot(d2,main="分解(by stl)")

par(opar)

cat("通过decompose函数得到季节项:",d1$seasonal[1:12],"\n")
## 通过decompose函数得到季节项: -16.56681 -12.22014 -5.965972 1.506528 7.281528 12.18236 13.83153 12.38903 7.375694 0.9131944 -7.470972 -13.25597
cat("通过decompose函数得到趋势项:",d1$trend,"\n")
## 通过decompose函数得到趋势项: NA NA NA NA NA NA 13.27083 13.10417 12.9375 12.85833 12.91667 13.0375 13.06667 12.9875 12.99583 12.99583 12.77917 12.6875 12.675 12.67917 12.85417 12.96667 12.90833 12.80833 12.8875 13.1125 13.1375 13.1 13.2 13.15 13.04583 13.0875 13.0875 13.0625 13.07917 13.03333 12.92083 12.7875 12.875 13.05833 13.03333 13.0375 13.19583 13.28333 13.15833 13.01667 12.975 13.03333 13.175 13.2625 13.22917 13.1 13.10417 13.15833 12.93333 12.57917 12.5625 12.70833 12.75417 12.84583 12.9625 13.02917 13.07083 13.09167 12.95417 12.83333 NA NA NA NA NA NA

指数平滑法拟合

#指数平滑
temfit.Holt <- HoltWinters(x.tem)
temfit.Holt
## Holt-Winters exponential smoothing with trend and additive seasonal component.
## 
## Call:
## HoltWinters(x = x.tem)
## 
## Smoothing parameters:
##  alpha: 0.004564983
##  beta : 0.754434
##  gamma: 0.4357686
## 
## Coefficients:
##              [,1]
## a    12.677240370
## b     0.003852541
## s1  -16.888229716
## s2  -12.352126979
## s3   -5.557169525
## s4    1.940114696
## s5    7.588304488
## s6   13.048077291
## s7   15.652727482
## s8   12.923083390
## s9    8.530897478
## s10   0.584495494
## s11  -8.289664171
## s12 -13.170275433
cat("当前点的水平、趋势、季节部分(系数)分别为",temfit.Holt$alpha,temfit.Holt$beta,temfit.Holt$gamma)
## 当前点的水平、趋势、季节部分(系数)分别为 0.004564983 0.754434 0.4357686
#ets
temfit.ets=ets(x.tem,model="AAA")
temfit.ets
## ETS(A,Ad,A) 
## 
## Call:
##  ets(y = x.tem, model = "AAA") 
## 
##   Smoothing parameters:
##     alpha = 1e-04 
##     beta  = 1e-04 
##     gamma = 2e-04 
##     phi   = 0.9515 
## 
##   Initial states:
##     l = 13.7512 
##     b = -0.0483 
##     s = -13.2145 -7.6068 0.9827 7.3885 12.6392 13.969
##            12.1381 7.3271 1.1773 -5.9342 -12.2006 -16.6658
## 
##   sigma:  1.3955
## 
##      AIC     AICc      BIC 
## 372.5193 385.4250 413.4993

从估计值中可以看出，alpha的值相对较低，表明当前时间点的水平估计是基于最近的观察和更远的过去的一些观察。斜率 beta 的估计值较高，说明其基于最近的观察。gamma的值较低，表明对季节性的估计基于最近的观察和更远的过去的一些观察。

  1. 预测未来12个月的月平均气温
#(3)
#指数平滑预测
plot(temfit.Holt)

temforecast.Holt<-forecast(temfit.Holt,h=12)
plot(temforecast.Holt)

cat("预测未来12个月的月平均气温(预测均值):",temforecast.Holt$mean)
## 预测未来12个月的月平均气温(预测均值): -4.207137 0.3328185 7.131628 14.63277 20.28481 25.74843 28.35694 25.63114 21.24281 13.30026 4.429954 -0.4468046

第四题

  1. 某城市1980年1月至1995年8月每月屠宰生猪数量如data2.txt所示
  1. 绘制时序图，直观考察序列的确定性特征
#(1)
data2 <- as.matrix(read.table("data2.txt",header=FALSE,sep="\t"))
data2 <- as.vector(t(data2))
data2 <- data2[-(which(is.na(data2)==TRUE))]
x.pig<-ts(data2,start=c(1980,1),frequency = 12)
plot(x.pig,ylab="Monthly pigs",main="每月屠宰生猪数量")

monthplot(x.pig, xlab="", ylab="")

seasonplot(x.pig, year.labels="TRUE", main="")

从时序图来看，数据具有明显的季节性和趋势性。不同年的同一个月宰猪数量呈先上升后下降再上升的趋势，但均值相差不大；总体的宰猪数量经历了波动平稳（1980-1985）-下降趋于波动平稳（1985-1990）-上升趋于波动平稳（1990-1995），从月度图和季节图也可以看出这点。

  1. 选择合适的方法对序列进行拟合（滑动平均，指数平滑等）
#(2)
#平滑处理
opar <- par(no.readonly = TRUE)
par(mfrow=c(2,2))
ylim <- c(min(x.pig)-100,max(x.pig)+100)
# 趋势图
plot(x.pig,main="Raw time series")
#ma
plot(ma(x.pig,10),main="Simple Moving Averages (k=10)",ylab="monthly pig",ylim=ylim)
#SMA
plot(SMA(x.pig,n=15),ylab="monthly pig",main="Simple Moving Averages (k=15)")
#rollmean
plot(rollmean(x.pig,20),ylab="monthly pig",main="Simple Moving Averages (k=20)",ylim=ylim)

par(opar)

#decompose和stl函数拟合趋势性和季节性(滑动平均去趋势)
par(mfrow=c(1,2))
d1=decompose(x.pig)
plot(d1)

d2=stl(x.pig,'periodic')
plot(d2,main="分解(by stl)")

par(opar)

cat("通过decompose函数得到季节项:",d1$seasonal[1:12],"\n")
## 通过decompose函数得到季节项: -10602.19 -7338.939 1882.524 -2662.636 2967.212 1665.105 5844.116 2155.283 -1546.647 991.9275 -181.9252 6826.175
cat("通过decompose函数得到趋势项:",d1$trend,"\n")
## 通过decompose函数得到趋势项: NA NA NA NA NA NA 89993.71 90404.33 93200.75 95599.5 95493.33 95583.96 95474.25 95038.42 94946.87 94519.04 94489.42 95271.54 95609.46 95796.33 96131.71 96177.96 95886.67 95599.5 95368.54 95113.42 95301 95037.62 94143.5 94135.17 94628.12 94916.33 95480.87 95860.71 96030 96778.54 97179.67 97549.29 97761.62 98322.46 99805.37 100182.8 100302.7 101151.6 101636.5 101772.6 102324.2 102937.5 103787.1 105094.4 105588.2 106024.6 106323.2 106482.2 107402.5 108088.6 108443.3 109115.2 109468.5 108728.6 107612.2 105640.5 103914.9 103206.5 102201 100339 98202.75 96360.96 93846.88 90979.08 88191.46 85876.38 84375.96 83356.62 82237.71 80477.79 78386.12 76833.25 75527.62 73879.92 73186.92 73272.25 72334.92 71095.17 70670.08 70367.96 70526.67 71344.87 72498.17 73626.96 74263.29 75344.67 76595.42 76768.71 77210.54 78043.12 77807.88 78083.21 78594.33 77947.75 77302.96 76809.42 76350.83 75993.33 75413.33 75240.67 75805.71 76295.71 76362.54 76177.38 76012.75 76136.37 76282.83 76348.38 76799.87 77383.67 77460.83 77917.5 79148.71 80380.12 81687.38 82540.75 82611.75 83197.25 83943.58 84488.33 85116.17 85334.12 85328.87 85354.96 85288.08 84786.79 84799.17 85665.96 86823.83 87911.75 89671.5 91514.5 92522.25 93761.83 95420.92 97162.38 98250.21 99752.04 100694.2 100232.8 99989.88 99428.5 97625.25 96265.33 95330.08 93796.67 92879.5 92180.62 91201.38 89957.83 89610.04 90511.88 91611.79 92728.04 93689.83 94234.75 94741.71 95313.04 96113.17 96870.42 97531.33 98541.62 98445.79 97664.96 97727.42 97969.42 98122.37 98392.04 98771.58 99308.83 99480.25 98876.92 98888.33 99563.87 99556.54 99453.04 NA NA NA NA NA NA
#指数平滑
pigfit.Holt <- HoltWinters(x.pig)
pigfit.Holt
## Holt-Winters exponential smoothing with trend and additive seasonal component.
## 
## Call:
## HoltWinters(x = x.pig)
## 
## Smoothing parameters:
##  alpha: 0.3074109
##  beta : 0.003657678
##  gamma: 0.4375144
## 
## Coefficients:
##             [,1]
## a   100572.16590
## b      211.71046
## s1     486.50198
## s2   -2040.98232
## s3   -4953.90627
## s4    1931.00406
## s5  -14281.72507
## s6  -11711.73657
## s7    3014.95841
## s8   -7888.68734
## s9    4979.71026
## s10    -88.60042
## s11  -5717.59270
## s12   -254.48340
cat("当前点的水平、趋势、季节部分(系数)分别为",pigfit.Holt$alpha,pigfit.Holt$beta,pigfit.Holt$gamma)
## 当前点的水平、趋势、季节部分(系数)分别为 0.3074109 0.003657678 0.4375144
#ets
pigfit.ets=ets(x.pig,model="AAA")
pigfit.ets
## ETS(A,Ad,A) 
## 
## Call:
##  ets(y = x.pig, model = "AAA") 
## 
##   Smoothing parameters:
##     alpha = 0.3037 
##     beta  = 1e-04 
##     gamma = 1e-04 
##     phi   = 0.9386 
## 
##   Initial states:
##     l = 91441.9326 
##     b = 754.24 
##     s = 6826.26 -181.714 992.3283 -1546.933 2155.48 5844.166
##            1665.114 3222.642 -2663.646 1882.659 -7338.432 -10857.92
## 
##   sigma:  9325.233
## 
##      AIC     AICc      BIC 
## 4439.453 4443.500 4497.709

从估计值中可以看出，alpha的值较低，表明当前时间点的水平估计是基于最近的观察和更远的过去的一些观察。斜率 beta 的估计值趋近于0，基本不会随着时间序列更新，这一点从时序图中也能直观感受到，水平序列变化很大，但是趋势大致相同。说明其基于最近的观察。gamma的值较低，表明对季节性的估计基于最近的观察和更远的过去的一些观察。

  1. 进行1年的预测
#(3)
#指数平滑预测
plot(pigfit.Holt)

pigforecast.Holt<-forecast(pigfit.Holt,h=12)
plot(pigforecast.Holt)

cat("预测未来12个月的月平均气温(预测均值):",pigforecast.Holt$mean)
## 预测未来12个月的月平均气温(预测均值): 101270.4 98954.6 96253.39 103350 87348.99 90130.69 105069.1 94377.16 107457.3 102600.7 97183.39 102858.2

第五题

  1. 读入hefei.csv数据，画出时序图,采用lubridate和dplyr包对合肥平均温度数据，按课本的例子用不同的方法进行分解。
#读入数据，画时序图
data3 <- read.csv("hefei.csv", head=T, sep = ",")
data3$CST <- ymd(data3$CST) # 利用 lubridate 包中的 ymd 函数来转换时间
plot(data3$CST, data3$Mean.TemperatureC, xlab = "time", ylab = "Mean Temperature",main="合肥平均温度图")

#数据处理得到2007-2017月平均温度数据
y <- select(data3,CST,Mean.TemperatureC)
y1 <- group_by(y,year(CST),month(CST)) 
d<- summarise(y1,month_mean=mean(Mean.TemperatureC,na.rm=TRUE))   #外层按照年，内层按照月份分组
## `summarise()` has grouped output by 'year(CST)'. You can override using the
## `.groups` argument.
data3.group <- c(d$month_mean)
data3.group <- ts(data3.group,start = c(2007,1), frequency=12)
#发现2017年8月只有三天的数据，用这三天数据平均代替该月平均值

#for(i in 1:12){
  #try <- filter(y,year(CST)==2008,month(CST)==i)
  #print(mean(try$Mean.TemperatureC))
#}           #检验上面分组计算是否正确
#得到季节项，画图
#用相同的月得到十二个季节项
get.season <- function(yd,n_NA){ # input: Detrended series,n_NA为缺失值个数
  ymat <- matrix(c(yd,rep(NA,n_NA)), byrow=TRUE, ncol=12)
  #print(ymat)   #检查
  ## season
  seas0 <- apply(ymat, 2, mean, na.rm=TRUE)
  return(seas0)
}

## 画去掉趋势项后的序列、季节项和随机项,返回随机项
plot.season <- function(yd, seas0){ # yd: Detrended series 
  ## season
  n <- end(yd)[1]-start(yd)[1]  #n为年数
  begin_year<-start(yd)[1] #begin_year为起始年份
  seas <- rep(seas0, n+1)
  seas <- ts(seas, start=start(yd), frequency=12)
  seas <- seas[1:length(yd)]
  seas <- ts(seas, frequency=12, start=start(yd))
  
  ## Error
  r <- c(yd) - seas
  r <- ts(r, frequency=12, start=start(yd))

  plot(yd, type='l', lwd=2,main="Detrended Series(black), Seasonal(red) and Random(cyan)",xlim=c(begin_year,begin_year+n), ylab="")
  abline(v=(begin_year):(begin_year+n), col="gray")
  abline(h=0,col= "grey")
  lines(seas, type="l", col="red")   #季节项
  lines(r, type="l", col="cyan")    #随机项
  
  return(r)
} 
#年平均法分解
y.group.mean <- data3.group  #重记该时间数据
ymat.mean <- matrix(c(y.group.mean,rep(NA,4)), byrow=TRUE, ncol=12)

cols <- rainbow(20)

#算趋势项
tr0.mean <- apply(ymat.mean, 1, mean, na.rm=TRUE)  #趋势项
tr.mean <- rep(tr0.mean,each=12)
n_NA <- length(tr.mean)-length(y.group.mean)   #缺失值个数
tr.mean <- tr.mean[1:length(y.group.mean)]   
tr.mean <- ts(tr.mean, frequency=12, start=start(y.group.mean))
cat("按照年平均得到趋势项为:",tr0.mean,"\n")
## 按照年平均得到趋势项为: 17.39055 16.14113 16.35157 16.4811 16.17697 16.32266 16.7845 16.41683 16.40289 16.7574 18.45437
#算季节项
y.detrended.mean <- y.group.mean - tr.mean
seas0.mean <- get.season(y.detrended.mean,n_NA)
cat("按照年平均法得到季节项为:",seas0.mean,"\n")
## 按照年平均法得到季节项为: -13.51929 -10.9961 -5.900526 -0.05575556 5.501234 8.586669 12.10827 11.42792 6.974106 1.887116 -5.419228 -11.30966
y.r.mean <- plot.season(y.detrended.mean, seas0.mean)

cat("按照年平均法得到随机项为:",y.r.mean,"\n")
## 按照年平均法得到随机项为: -0.5486782 2.534116 0.2196503 -0.7681308 1.172729 -0.6438883 -0.9181796 0.4395916 -0.1646587 -0.4712178 -0.5713254 0.4352338 -2.105709 -2.627793 1.82391 -0.05204323 1.131828 -0.8611341 0.6860799 -0.3432456 1.351429 0.745945 0.3780955 0.5878804 -0.5742158 1.573094 -0.6768551 0.5041811 -0.4334533 1.728423 0.1530569 -0.06981698 0.5409865 1.954857 -3.26568 -0.7193361 0.7801301 0.6935692 -1.64509 -2.825344 -0.3049139 0.9322317 0.2815963 1.123239 0.3114614 -1.14241 1.004795 1.505977 -2.367355 0.3191261 -0.6635423 1.378784 0.3540531 0.1696933 0.06959558 -0.3468267 -0.151077 -0.2898942 3.142256 -0.8995717 0.2934074 -1.878288 -1.099554 1.933095 0.6277182 1.557337 1.891648 0.668774 -0.7634334 0.2095774 -0.9701 -1.754939 -0.5232685 -0.2884008 0.9224793 -0.262076 -0.1889577 -0.2378336 1.171746 2.465001 -0.6252706 -0.6716147 0.2347294 -1.281292 2.844403 -1.456444 1.677248 0.2055956 0.6303267 0.1965046 -1.234776 -2.231843 -0.7575991 0.5347665 0.7690676 -0.4620077 2.277692 0.8074901 0.4653748 -0.4471374 0.06361527 -0.3895617 -2.059552 -0.5082322 -0.1436654 0.3228938 -0.4836654 0.8099906 -0.238106 1.031801 0.9818353 0.865022 -2.194118 -0.9774022 -0.1882527 0.6856476 0.4018274 -1.192904 -0.2381726 1.778064 0.1617001 -0.7082709 -2.005455 -0.531946 -0.8588278 -1.47437 0.1470373 -1.882288
#线性回归直线拟合法
cols <- rainbow(20)

y.group.lm <- data3.group  #重记该时间数据
tt <- seq(length(y.group.lm))
lmr <- lm(y.group.lm ~ tt)

#趋势项
tr.trend.lm <- ts(predict(lmr,newdata=list(tt=seq(length(y.group.lm)))),frequency=12, start=start(y.group.lm))   
cat("按照线性回归拟合得到趋势项为:",tr.trend.lm,"\n")
## 按照线性回归拟合得到趋势项为: 15.97796 15.98844 15.99891 16.00939 16.01987 16.03035 16.04082 16.0513 16.06178 16.07226 16.08274 16.09321 16.10369 16.11417 16.12465 16.13513 16.1456 16.15608 16.16656 16.17704 16.18751 16.19799 16.20847 16.21895 16.22943 16.2399 16.25038 16.26086 16.27134 16.28181 16.29229 16.30277 16.31325 16.32373 16.3342 16.34468 16.35516 16.36564 16.37611 16.38659 16.39707 16.40755 16.41803 16.4285 16.43898 16.44946 16.45994 16.47042 16.48089 16.49137 16.50185 16.51233 16.5228 16.53328 16.54376 16.55424 16.56472 16.57519 16.58567 16.59615 16.60663 16.6171 16.62758 16.63806 16.64854 16.65902 16.66949 16.67997 16.69045 16.70093 16.7114 16.72188 16.73236 16.74284 16.75332 16.76379 16.77427 16.78475 16.79523 16.80571 16.81618 16.82666 16.83714 16.84762 16.85809 16.86857 16.87905 16.88953 16.90001 16.91048 16.92096 16.93144 16.94192 16.95239 16.96287 16.97335 16.98383 16.99431 17.00478 17.01526 17.02574 17.03622 17.04669 17.05717 17.06765 17.07813 17.08861 17.09908 17.10956 17.12004 17.13052 17.14099 17.15147 17.16195 17.17243 17.18291 17.19338 17.20386 17.21434 17.22482 17.2353 17.24577 17.25625 17.26673 17.27721 17.28768 17.29816 17.30864
#季节项
y.detrended.lm <- y.group.lm - tr.trend.lm[1:length(y.group.lm)] #去掉趋势项后的序列
seas0.lm <- get.season(y.detrended.lm,n_NA)   #季节项
cat("按照线性回归拟合得到季节项为:",seas0.lm,"\n")
## 按照线性回归拟合得到季节项为: -13.42774 -10.91502 -5.829928 0.004364232 5.550876 8.625833 12.13696 11.44613 6.869084 1.771617 -5.545204 -11.44611
y.r.lm <- plot.season(y.detrended.lm,seas0.lm)

cat("按照线性回归拟合得到随机项为:",y.r.lm,"\n")
## 按照线性回归拟合得到随机项为: 0.7723637 3.855158 1.540692 0.5529111 2.493771 0.6771536 0.4028623 1.760634 1.269135 0.962576 0.8624685 1.869028 -2.159822 -2.681906 1.769797 -0.1061559 1.077715 -0.9152468 0.6319672 -0.3973583 1.410068 0.8045842 0.4367348 0.6465197 -0.5436199 1.60369 -0.6462592 0.534777 -0.4028574 1.759019 0.1836528 -0.03922106 0.6843344 2.098205 -3.122332 -0.5759882 0.8145173 0.7279565 -1.610703 -2.790957 -0.2705266 0.966619 0.3159836 1.157626 0.4586006 -0.9952704 1.151934 1.653117 -2.762829 -0.07634866 -1.059017 0.9833095 -0.04142168 -0.2257814 -0.3258792 -0.7423014 -0.4337998 -0.572617 2.859534 -1.182294 -0.08211144 -2.253807 -1.475073 1.557576 0.2521994 1.181818 1.516129 0.2932551 -1.0262 -0.05318945 -1.232867 -2.017706 -0.5626839 -0.3278161 0.8830639 -0.3014913 -0.228373 -0.2772489 1.132331 2.425586 -0.5519339 -0.598278 0.3080661 -1.207955 2.311582 -1.989264 1.144427 -0.3272251 0.09750612 -0.336316 -1.767597 -2.764664 -1.177668 0.1146979 0.348999 -0.8820763 1.605203 0.1350021 -0.2071133 -1.119625 -0.6088728 -1.06205 -2.73204 -1.18072 -0.7034014 -0.2368423 -1.043401 0.2502545 -0.6818206 0.5880861 0.5381207 0.4213074 -2.637832 -1.421117 -0.6319672 0.2419331 0.07086482 -1.523866 -0.5691352 1.447101 1.28922 0.4192489 -0.8779356 0.5955737 0.268692 -0.3468505 1.274557 -0.7547684
#画趋势项与原序列的图
seas.lm <- ts(rep(seas0.lm, end(y.detrended.lm)[1]-start(y.detrended.lm)[1]+1),
                start=start(y.detrended.lm), frequency=12)
y.pred.lm <- tr.trend.lm + seas.lm   #拟合值
plot(y.group.lm, main="Series(black), Linear trend(red)",
     lwd=2,type="l", col="black")
lines(tr.trend.lm, col="red", type="l")    #趋势项的图

#二次曲线拟合法
cols <- rainbow(20)

y.group.qlm <- data3.group  #重记该时间数据
tt <- seq(length(y.group.qlm))
lmr <- lm(y.group.lm ~ tt+I(tt^2))

#趋势项
tr.trend.qlm <- ts(predict(lmr,newdata=list(tt=seq(length(y.group.qlm)))),frequency=12, start=start(y.group.qlm))   
cat("按照二次曲线拟合得到趋势项为:",tr.trend.qlm,"\n")
## 按照二次曲线拟合得到趋势项为: 16.81718 16.78801 16.75947 16.73156 16.70428 16.67762 16.6516 16.6262 16.60144 16.5773 16.5538 16.53092 16.50867 16.48705 16.46606 16.4457 16.42597 16.40687 16.3884 16.37056 16.35335 16.33676 16.32081 16.30548 16.29079 16.27672 16.26328 16.25047 16.2383 16.22675 16.21583 16.20554 16.19588 16.18684 16.17844 16.17067 16.16353 16.15701 16.15113 16.14587 16.14124 16.13725 16.13388 16.13114 16.12903 16.12755 16.1267 16.12648 16.12689 16.12793 16.12959 16.13189 16.13482 16.13837 16.14256 16.14737 16.15281 16.15888 16.16559 16.17292 16.18088 16.18947 16.19869 16.20854 16.21901 16.23012 16.24186 16.25422 16.26722 16.28084 16.2951 16.30998 16.32549 16.34163 16.35841 16.37581 16.39384 16.41249 16.43178 16.4517 16.47225 16.49343 16.51523 16.53767 16.56073 16.58442 16.60875 16.6337 16.65928 16.68549 16.71233 16.7398 16.7679 16.79663 16.82599 16.85598 16.88659 16.91784 16.94972 16.98222 17.01535 17.04912 17.08351 17.11853 17.15418 17.19047 17.22738 17.26491 17.30308 17.34188 17.38131 17.42137 17.46205 17.50337 17.54531 17.58789 17.63109 17.67492 17.71938 17.76448 17.8102 17.85655 17.90353 17.95114 17.99937 18.04824 18.09774 18.14786
#季节项
y.detrended.qlm <- y.group.qlm - tr.trend.qlm[1:length(y.group.qlm)] #去掉趋势项后的序列
seas0.qlm <- get.season(y.detrended.qlm,n_NA)   #季节项
cat("按照二次曲线拟合得到季节项为:",seas0.qlm,"\n")
## 按照二次曲线拟合得到季节项为: -13.45512 -10.94051 -5.854158 -0.01923601 5.527275 8.601603 12.11147 11.41875 6.924152 1.827314 -5.489508 -11.39105
y.r.qlm <- plot.season(y.detrended.qlm,seas0.qlm)

cat("按照二次曲线拟合得到随机项为:",y.r.qlm,"\n")
## 按照二次曲线拟合得到随机项为: -0.03948468 3.08107 0.8043646 -0.1456561 1.832965 0.05410711 -0.1824238 1.213108 0.674409 0.4018342 0.3357111 1.376255 -2.537426 -3.029302 1.45261 -0.3931349 0.8209444 -1.141809 0.4356132 -0.563504 1.18917 0.6101182 0.268701 0.5049182 -0.5776042 1.592362 -0.6349311 0.5687613 -0.3462168 1.838316 0.2856059 0.08538823 0.746639 2.17939 -3.022267 -0.457043 1.033528 0.9620709 -1.361484 -2.526634 0.008900247 1.26115 0.6256188 1.482365 0.7134833 -0.7290596 1.429473 1.941984 -2.381449 0.3125834 -0.6625329 1.387346 0.3701666 0.1933589 0.1008133 -0.3080569 -0.0769641 -0.2120053 3.223921 -0.8141306 0.3710133 -1.800682 -1.021949 2.0107 0.7053241 1.634943 1.969254 0.7463798 -0.6580364 0.3111983 -0.8722551 -1.66087 -0.1284394 0.09887628 1.302204 0.1100969 0.1756631 0.1192352 1.521263 2.806966 -0.263067 -0.3207392 0.5742768 -0.9530728 2.636322 -1.679629 1.438958 -0.04779817 0.3618289 -0.08709739 -1.533482 -2.545654 -1.058722 0.2147629 0.4301838 -0.8197716 1.729813 0.2369551 -0.1278164 -1.062985 -0.5748884 -1.050722 -2.743368 -1.214705 -0.8450029 -0.404876 -1.237867 0.02935619 -0.8479663 0.391732 0.3115584 0.1645368 -2.924811 -1.738304 -0.9793628 -0.1356709 -0.4219083 -2.050624 -1.129877 0.8523752 0.7416942 -0.1660372 -1.500982 -0.06523314 -0.4298752 -1.083178 0.5004691 -1.566617
#画趋势项与原序列的图
seas.qlm <- ts(rep(seas0.qlm, end(y.detrended.qlm)[1]-start(y.detrended.qlm)[1]+1),
                start=start(y.detrended.qlm), frequency=12)
y.pred.qlm <- tr.trend.qlm + seas.qlm   #拟合值
plot(y.group.qlm, main="Series(black), Quadratic trend(red)",
     lwd=2,type="l", col="black")
lines(tr.trend.qlm, col="red", type="l")    #趋势项的图