非平稳时间序列趋势分析

有些时间序列具有非常显著的趋势，有时我们分析的目的就是要找到序列中的这种趋势，并利用这种趋势对未来的发展作出合理的预测。

趋势拟合法

趋势拟合法就是把时间作为自变量，相应的序列观测值做为因变量，建立序列值随时间变化的回归模型的方法。根据序列所表现出的线性或者非线性特征，拟合方法又可以具体分为线性拟合曲线拟合。

线性拟合

如果长期趋势呈现出线性特征，那我们可以用线性模型来拟合。

R语言使用lm函数拟合线性趋势。

lm(Y ~ a + X1 + X2 +…+Xn,data = )

-Y:响应变量

-a:指定是否需要常数项

（1）a = 1，模型有非零常数项，这是默认设置。

（2）a = 0，模型不需要常数项。

-X1,…,Xn:自变量。

-data:数据框名。如果自变量和响应变量不是独立输入变量而是共同存在于某个数据框中，则需要指定数据框名。

#读入数据
x<-c(8444,9215,8879,8990,8115,9457,8590,9294,8997,9574,9051,9724,9120,
+ 10143,9746,10074,9578,10817,10116,10779,9901,11266,10686,10961,10121,
+ 11333,10677,11325,10698,11624,11502,11393,10609,12077,11376,11777,
+ 11225,12231,11884,12109)
#构造时间变量
t<-c(1:40)
#拟合回归模型
x.fit<-lm(x~t)
#查看拟合信息
summary(x.fit)
Call:
lm(formula = x ~ t)

Residuals:
Min      1Q  Median      3Q     Max
-853.06 -329.64   63.54  314.46  794.04

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 8491.765    137.964   61.55   <2e-16 ***
t             90.009      5.864   15.35   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 428.1 on 38 degrees of freedom
Multiple R-squared:  0.8611,    Adjusted R-squared:  0.8575
F-statistic: 235.6 on 1 and 38 DF,  p-value: < 2.2e-16

#绘制拟合效果图
x<-ts(x)
plot(x)

abline(lm(x~t),col=2)

曲线拟合

如果长期趋势呈现出非线性特征，那么我们可以用曲线典线模型来拟合它。对曲线模型进行参数估计时，能转换成线性模型的都转化成线性模型，用线性最小二乘法进行参数估计，实在不能转化成线性模型的，就用迭代法进行参数估计。





R语言针对非线性趋势的拟合也分为两类：一类可以写成关于时间t的多项式，这时仍然可以用lm函数拟合。另一类无法通过适当的变换变成线性回归模型，只能通过非线性回归解决，这时要用nls函数。

nls(Y ~ f(x1,…,xn),data=,start=)

-Y:响应变量。

-X1,…,Xn：自变量

-f:非线性函数

-data:数据框名

-start:如果需要利用迭代计算未知参数，可以指定迭代初始值。

1949-2008年化肥产量序列进行曲线拟合’

lm函数拟合

a<-read.table("D:/R-TT/book4/4R/data/file12.csv",sep=",",header = T)
x<-ts(a$output,start=1949)
#lm函数拟合
t1<-c(1:60)
t2<-t1^2
x.fit1<-lm(x~t1+t2)
summary(x.fit1)
Call:
lm(formula = x ~ t1 + t2)

Residuals:
Min      1Q  Median      3Q     Max
-532.12 -164.92   24.68  105.51  716.37

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 319.0255   105.3371   3.029  0.00369 **
t1          -57.7690     7.9679  -7.250 1.22e-09 ***
t2            2.3551     0.1266  18.601  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 263 on 57 degrees of freedom
Multiple R-squared:  0.9755,    Adjusted R-squared:  0.9746
F-statistic:  1133 on 2 and 57 DF,  p-value: < 2.2e-16

nls函数拟合

x.fit2<-nls(x~a+b*t1+c*t1^2,start = list(a=1,b=1,c=1))
> summary(x.fit2)

Formula: x ~ a + b * t1 + c * t1^2

Parameters:
Estimate Std. Error t value Pr(>|t|)
a 319.0255   105.3371   3.029  0.00369 **
b -57.7690     7.9679  -7.250 1.22e-09 ***
c   2.3551     0.1266  18.601  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 263 on 57 degrees of freedom

Number of iterations to convergence: 1
Achieved convergence tolerance: 3.059e-08

根据输出结果可以知道，这两个函数得到的拟合结果完全一致，1949-2008年我国化肥序列的拟合模型为：

Xt = 319.0255 - 57.769t +2.3551t^2 N（0，263^2）

y<-predict(x.fit2)
y<-ts(y,start = 1949)
plot(x,type = "p")
lines(y,col=2,lwd=2)

平滑法

平滑法是进行趋势分析和预测时常用的一种方法，它是利用修匀技术，消弱短期随机波动的序列的影响，使序列平滑化，从而显示出变化的规律。它具有调节灵活、计算简便的特征，广泛应用于计量经济、人口研究等诸多领域。根据所有的平滑技术的不同，平滑法又可以具体分为移动平均法与指数平滑法。

移动平均法

移动平均法的基本思想是对于一个时间序列{Xt},我们可以假定在一个比较短的时间间隔内，序列的预测值是比较平稳的，它们之间的差异主要是由随机波动造成的，根据这种假定，我们可以使用一定时间间隔的平均值作为下一期的估计值。

移动平均的期数对原序列的修匀效果影响很大，要确定移动平均的期数，一般会从如下三个方面加以考虑：

（1）事件的发展有无周期性。

（2）对趋势平滑性的要求。

（3）对趋势反映近期灵敏度的要求。

在R语言中，TTR程序包中的SMA函数是专门用来作简单移动平均趋势拟合的函数。

SMA(x,n)

-x:需要做简单的移动平均的序列名

-n:移动平均期数

对北京市1949年-1998年每年最高气温序列进行进行5期移动平均拟合。

library(TTR)
a<-read.table("D:/R-TT/book4/4R/data/file6.csv",",",header = T)
x<-ts(a$temp,start = 1949)
x.ma<-SMA(x,n=5)
plot(x,type = "o")
lines(x.ma,col=2,lwd=2)

指数平滑法

移动平均法实际上就是用一个简单的加权平均数作为某一期的趋势的估计值。

原理为近期的结果对现在的影响大些，远期的结果对现在的影响小些。

各期权重随时间间隔的增大而呈指数衰减，这就是指数平滑法的基本思想。

1.简单的指数平滑



2.Holt两参数指数平滑

适用于对含有线性趋势的序列进行修匀。它的基本思想是假定序列有一个比较固定的线性趋势——每期都递增r或者递减r,那么第t期的估计值就应该等于第t-1期的观察值加上每期固定的趋势变动。





3.Holt-Winters三参数指数平滑

假定要进行指数平滑的序列为{Xt},{Xt}序列既含有趋势又含有季节。这个季节因子可以随每年的具体情况波动，因为引入了季节因子，所以构建了三参数模型。

在R语言中，HoltWinters函数可以完成上述三种平滑趋势拟合。

HoltWinters(x,alpha=,beta=,gamma=,seasonal=)

-x:要进行指数平滑的序列名

-alpha:随机波动部分的参数

-beta:趋势部分的参数

-gamma:季节部分的参数

这三个指数联合起来，确定要拟合的指数平滑模型类型：

（1）当alpha不指定时，beta = F,gamma=F,表示拟合简单指数平滑模型。

（2）当alpha和beta不指定时，gamma=F时，表示拟合Holt两参数指数平滑模型。

（3）当三个参数都不指定时，表示拟合Holt-Winters三参数指数平滑模型。

-seasonal:当既含有季节的又含有趋势时，指定季节与趋势的关系。

seasonal = “additive”表示加法关系， 这是系统默认选项。seasonal=”multiplicative”表示乘法关系。

对1964-1999年中国纱年产量序列分析,并预测未来10年的

#读入序列
a<-read.table("D:/R-TT/book4/4R/data/file4.csv",sep=",",header = T)
x<-ts(a$output,start = 1964)
#进行Holt两参数平滑
x.fit<-HoltWinters(x,gamma = F)
x.fit
Holt-Winters exponential smoothing with trend and without seasonal component.

Call:
HoltWinters(x = x, gamma = F)

Smoothing parameters:
alpha: 0.855644
beta : 0.158537
gamma: FALSE

Coefficients:
[,1]
a 565.55301
b  12.29066

#绘制Holt两参数指数平滑拟合效果图
plot(x.fit)

#预测序列并绘制预测效果图
library(forecast)
x.fore<-forecast(x.fit,h=10)
x.fore
Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
2000       577.8437 545.8643 609.8230 528.9355 626.7519
2001
c512
590.1343 545.1051 635.1636 521.2681 659.0006
2002       602.4250 544.7497 660.1003 514.2182 690.6318
2003       614.7157 544.3115 685.1198 507.0417 722.3896
2004       627.0063 543.6025 710.4101 499.4512 754.5614
2005       639.2970 542.5375 736.0565 491.3161 787.2779
2006       651.5876 541.0751 762.1002 482.5733 820.6020
2007       663.8783 539.1960 788.5606 473.1931 854.5635
2008       676.1690 536.8922 815.4458 463.1635 889.1744
2009       688.4596 534.1622 842.7570 452.4821 924.4371
plot(x.fore)

对1962年1月和1975年12月平均每头奶牛月产奶量序列进行Holt-Winters三参数指数平滑

#读入序列
b<-read.table("D:/R-TT/book4/4R/data/file5.csv",sep=",",header = T)
x<-ts(b$milk,start = c(1962,1),frequency = 12)
#进行Holt-winters三参数指数平滑
x.fit<-HoltWinters(x)
x.fit
Holt-Winters exponential smoothing with trend and additive seasonal component.

Call:
HoltWinters(x = x)

Smoothing parameters:
alpha: 0.68933
beta : 0
gamma: 0.8362592

Coefficients:
[,1]
a   885.775547
b     1.278118
s1  -16.743296
s2  -59.730034
s3   47.492731
s4   56.203890
s5  115.537545
s6   84.554817
s7   39.580306
s8   -4.702033
s9  -54.554684
s10 -51.582594
s11 -85.953466
s12 -42.907363

#绘制Holt-winters三参数指数平滑拟合效果图
plot(x.fit)

#预测序列并绘制预测效果图
x.fore<-forecast(x.fit,h=24)
plot(x.fore)

月度产奶量Holt-Winters三参数指数平滑序列预测图