Time Series Analysis Chapter 4
1 非平稳时间序列的确定性分析
1.1 确定性因素分解
1.1.1 趋势分析
1.1.1.1 趋势拟合法
1.1.1.1.1 线性拟合：线性回归
xt=a+bt+It,t=1,2,...,40E(It)=0,Var(It)=σ2
# 这是澳洲1981年-1990年每季度的消费支出时间序列
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由输出效果可知，线性拟合模型为
xt=8491.765+90.009t+ϵt,ϵt∼WN(0,428.12)
# 绘制拟合效果图
x<-ts(x)
plot(x)
abline(lm(x~t),col=2)
1.1.1.1.2 曲线拟合：非线性回归
能转换成线性模型的，均转换成线性模型，用最小二乘法进行参数估计；不能转换成线性模型的，运用迭代法进行参数估计。
xt=a0+a1t+a2t2,t=1,2,...,60
# 这是我国1949年-2008年化肥产量序列
a<-read.table("file12.csv",sep=",",header = T)
x<-ts(a$output,start=1949)
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函数进行拟合，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由输出结果，这两个函数的拟合结果完全一致，故建立模型为
xt=319.0255−57.7690t+2.3551t2+ϵt,ϵt∼WN(0,2632)
# 绘图
y<-predict(x.fit2) # 拟合值赋值给y
y<-ts(y,start = 1949)
plot(x,type = "p")
lines(y,col=2,lwd=2)
1.1.1.2 平滑法
1.1.1.2.1 移动平均法
# 这是北京市1949年-1998年每年最高气温序列，对其进行 5 期移动平均拟合
library(TTR)
a<-read.table("file6.csv",",",header = T)
x<-ts(a$temp,start = 1949)
x.ma<-SMA(x,n=5) # 作 5 期移动平均拟合
plot(x,type = "o")
lines(x.ma,col=2,lwd=2)
1.1.1.2.2 指数平滑法
# 这是1964年-1999年我国纱年产量序列，对其进行 Holt 两参数指数平滑，并预测未来 10 年的中国纱产量。
a<-read.table("file4.csv",sep=",",header = T)
x<-ts(a$output,start = 1964)
x.fit<-HoltWinters(x,gamma = F) # 进行Holt两参数平滑，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对于 HoltWinters() 函数，有如下参数搭配：
若 alpha 不指定，beta=F、gamma=F 时，表示拟合简单指数平滑模型。
若 alpha、beta 不指定，gamma=F 时。表示拟合 Holt 两参数指数平滑模型。
若三个指数均不指定时，表示拟合 Holt-Winters 三参数指数平滑模型。
seasonal 表示既含有季节又含有趋势时，指定季节与趋势的关系。’seasonal=“addictive”’默认参数，表示加法关系；’seasonal=“multiplicative”’表示乘法关系。
plot(x.fit) # 绘制Holt两参数指数平滑拟合效果图
# 预测序列
library(forecast)## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoox.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 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)
# 这是1961年1月-1975年12月平均每头奶牛月产奶量序列，对其进行 Holt-Winters 三参数指数平滑，并预测未来 2 年平均每头奶牛的月度产奶量。
b<-read.table("file5.csv",sep=",",header = T)
x<-ts(b$milk,start = c(1962,1),frequency = 12)
x.fit<-HoltWinters(x) # 进行Holt-winters三参数指数平滑
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.907363plot(x.fit) # 绘制Holt-winters三参数指数平滑拟合效果图
#预测序列
x.fore<-forecast(x.fit,h=24)
x.fore## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## Jan 1976 870.3104 859.7804 880.8403 854.2062 886.4145
## Feb 1976 828.6017 815.8124 841.3911 809.0421 848.1614
## Mar 1976 937.1026 922.3970 951.8082 914.6124 959.5929
## Apr 1976 947.0919 930.6925 963.4914 922.0111 972.1727
## May 1976 1007.7037 989.7697 1025.6377 980.2760 1035.1314
## Jun 1976 977.9991 958.6518 997.3463 948.4100 1007.5882
## Jul 1976 934.3027 913.6386 954.9668 902.6997 965.9057
## Aug 1976 891.2985 869.3966 913.2003 857.8025 924.7945
## Sep 1976 842.7239 819.6506 865.7973 807.4363 878.0115
## Oct 1976 846.9741 822.7860 871.1623 809.9816 883.9667
## Nov 1976 813.8814 788.6276 839.1352 775.2590 852.5037
## Dec 1976 858.2056 831.9294 884.4818 818.0195 898.3916
## Jan 1977 885.6478 857.5350 913.7605 842.6530 928.6425
## Feb 1977 843.9392 814.9045 872.9739 799.5344 888.3439
## Mar 1977 952.4400 922.5118 982.3683 906.6687 998.2114
## Apr 1977 962.4293 931.6334 993.2252 915.3310 1009.5276
## May 1977 1023.0411 991.4013 1054.6809 974.6522 1071.4300
## Jun 1977 993.3365 960.8748 1025.7982 943.6906 1042.9824
## Jul 1977 949.6401 916.3767 982.9035 898.7682 1000.5120
## Aug 1977 906.6359 872.5897 940.6820 854.5668 958.7049
## Sep 1977 858.0613 823.2500 892.8726 804.8221 911.3006
## Oct 1977 862.3115 826.7515 897.8715 807.9272 916.6959
## Nov 1977 829.2188 792.9255 865.5120 773.7130 884.7245
## Dec 1977 873.5430 836.5310 910.5550 816.9380 930.1480# 绘制预测效果图
plot(x.fore)
1.1.2 综合分析：确定性因素分解
# 这是1993年-2000年我国社会消费品零售总额序列，对其进行确定性因素分解
c<-read.table("file14.csv",sep=",",header = T)
x<-ts(c$sales,start = c(1993,1),frequency = 12)
plot(x)
# 确定性因素分解
x.fit<-decompose(x,type="mult")
x.fit## $x
## Jan Feb Mar Apr May Jun Jul Aug Sep Oct
## 1993 977.5 892.5 942.3 941.3 962.2 1005.7 963.8 959.8 1023.3 1051.1
## 1994 1192.2 1162.7 1167.5 1170.4 1213.7 1281.1 1251.5 1286.0 1396.2 1444.1
## 1995 1602.2 1491.5 1533.3 1548.7 1585.4 1639.7 1623.6 1637.1 1756.0 1818.0
## 1996 1909.1 1911.2 1860.1 1854.8 1898.3 1966.0 1888.7 1916.4 2083.5 2148.3
## 1997 2288.5 2213.5 2130.9 2100.5 2108.2 2164.7 2102.5 2104.4 2239.6 2348.0
## 1998 2549.5 2306.4 2279.7 2252.7 2265.2 2326.0 2286.1 2314.6 2443.1 2536.0
## 1999 2662.1 2538.4 2403.1 2356.8 2364.0 2428.8 2380.3 2410.9 2604.3 2743.9
## 2000 2774.7 2805.0 2627.0 2572.0 2637.0 2645.0 2597.0 2636.0 2854.0 3029.0
## Nov Dec
## 1993 1102.0 1415.5
## 1994 1553.8 1932.2
## 1995 1935.2 2389.5
## 1996 2290.1 2848.6
## 1997 2454.9 2881.7
## 1998 2652.2 3131.4
## 1999 2781.5 3405.7
## 2000 3108.0 3680.0
##
## $seasonal
## Jan Feb Mar Apr May Jun Jul
## 1993 1.0439030 0.9939439 0.9592626 0.9397644 0.9438897 0.9588798 0.9286603
## 1994 1.0439030 0.9939439 0.9592626 0.9397644 0.9438897 0.9588798 0.9286603
## 1995 1.0439030 0.9939439 0.9592626 0.9397644 0.9438897 0.9588798 0.9286603
## 1996 1.0439030 0.9939439 0.9592626 0.9397644 0.9438897 0.9588798 0.9286603
## 1997 1.0439030 0.9939439 0.9592626 0.9397644 0.9438897 0.9588798 0.9286603
## 1998 1.0439030 0.9939439 0.9592626 0.9397644 0.9438897 0.9588798 0.9286603
## 1999 1.0439030 0.9939439 0.9592626 0.9397644 0.9438897 0.9588798 0.9286603
## 2000 1.0439030 0.9939439 0.9592626 0.9397644 0.9438897 0.9588798 0.9286603
## Aug Sep Oct Nov Dec
## 1993 0.9260807 0.9814290 1.0074970 1.0472403 1.2694493
## 1994 0.9260807 0.9814290 1.0074970 1.0472403 1.2694493
## 1995 0.9260807 0.9814290 1.0074970 1.0472403 1.2694493
## 1996 0.9260807 0.9814290 1.0074970 1.0472403 1.2694493
## 1997 0.9260807 0.9814290 1.0074970 1.0472403 1.2694493
## 1998 0.9260807 0.9814290 1.0074970 1.0472403 1.2694493
## 1999 0.9260807 0.9814290 1.0074970 1.0472403 1.2694493
## 2000 0.9260807 0.9814290 1.0074970 1.0472403 1.2694493
##
## $trend
## Jan Feb Mar Apr May Jun Jul Aug
## 1993 NA NA NA NA NA NA 1028.696 1048.900
## 1994 1153.913 1179.492 1208.621 1240.533 1275.733 1316.087 1354.700 1385.483
## 1995 1537.554 1567.687 1597.308 1627.879 1659.350 1694.296 1726.138 1756.412
## 1996 1890.954 1913.637 1938.921 1966.329 1994.879 2028.796 2063.733 2092.137
## 1997 2190.733 2207.475 2221.812 2236.637 2251.825 2260.071 2272.325 2287.071
## 1998 2350.200 2366.608 2383.846 2400.158 2416.213 2434.837 2449.933 2464.292
## 1999 2513.642 2521.579 2532.308 2547.687 2561.737 2578.554 2594.675 2610.475
## 2000 2707.971 2726.379 2746.162 2768.446 2793.929 2818.963 NA NA
## Sep Oct Nov Dec
## 1993 1069.542 1088.471 1108.496 1130.450
## 1994 1414.425 1445.429 1476.679 1507.108
## 1995 1787.517 1813.887 1839.679 1866.312
## 1996 2116.017 2137.537 2156.521 2173.546
## 1997 2297.142 2309.683 2322.567 2335.829
## 1998 2479.100 2488.579 2497.033 2505.433
## 1999 2630.912 2649.208 2669.550 2689.933
## 2000 NA NA NA NA
##
## $random
## Jan Feb Mar Apr May Jun Jul
## 1993 NA NA NA NA NA NA 1.0088882
## 1994 0.9897286 0.9917699 1.0069996 1.0039379 1.0079296 1.0151591 0.9947887
## 1995 0.9982198 0.9571982 1.0006930 1.0123394 1.0122310 1.0092784 1.0128538
## 1996 0.9671360 1.0048115 1.0000891 1.0037414 1.0081543 1.0106040 0.9854906
## 1997 1.0006939 1.0088390 0.9998115 0.9993280 0.9918727 0.9988758 0.9963425
## 1998 1.0391783 0.9804972 0.9969239 0.9987216 0.9932308 0.9962666 1.0048103
## 1999 1.0145206 1.0128044 0.9892766 0.9843682 0.9776685 0.9823163 0.9878519
## 2000 0.9815488 1.0351058 0.9972322 0.9885893 0.9999390 0.9785257 NA
## Aug Sep Oct Nov Dec
## 1993 0.9880930 0.9748693 0.9584809 0.9492950 0.9863775
## 1994 1.0022841 1.0057935 0.9916461 1.0047607 1.0099323
## 1995 1.0064677 1.0009573 0.9948091 1.0044711 1.0085729
## 1996 0.9891157 1.0032647 0.9975563 1.0140385 1.0323984
## 1997 0.9935731 0.9933992 1.0090249 1.0092977 0.9718346
## 1998 1.0142266 1.0041262 1.0114724 1.0142280 0.9845558
## 1999 0.9972656 1.0086157 1.0280362 0.9949349 0.9973544
## 2000 NA NA NA NA NA
##
## $figure
## [1] 1.0439030 0.9939439 0.9592626 0.9397644 0.9438897 0.9588798 0.9286603
## [8] 0.9260807 0.9814290 1.0074970 1.0472403 1.2694493
##
## $type
## [1] "multiplicative"
##
## attr(,"class")
## [1] "decomposed.ts"plot(x.fit)
#查看季节指数
x.fit$figure## [1] 1.0439030 0.9939439 0.9592626 0.9397644 0.9438897 0.9588798 0.9286603
## [8] 0.9260807 0.9814290 1.0074970 1.0472403 1.2694493# 绘制季节指数图
plot(x.fit$figure,type = "o")
# 查看趋势拟合值
x.fit$trend## Jan Feb Mar Apr May Jun Jul Aug
## 1993 NA NA NA NA NA NA 1028.696 1048.900
## 1994 1153.913 1179.492 1208.621 1240.533 1275.733 1316.087 1354.700 1385.483
## 1995 1537.554 1567.687 1597.308 1627.879 1659.350 1694.296 1726.138 1756.412
## 1996 1890.954 1913.637 1938.921 1966.329 1994.879 2028.796 2063.733 2092.137
## 1997 2190.733 2207.475 2221.812 2236.637 2251.825 2260.071 2272.325 2287.071
## 1998 2350.200 2366.608 2383.846 2400.158 2416.213 2434.837 2449.933 2464.292
## 1999 2513.642 2521.579 2532.308 2547.687 2561.737 2578.554 2594.675 2610.475
## 2000 2707.971 2726.379 2746.162 2768.446 2793.929 2818.963 NA NA
## Sep Oct Nov Dec
## 1993 1069.542 1088.471 1108.496 1130.450
## 1994 1414.425 1445.429 1476.679 1507.108
## 1995 1787.517 1813.887 1839.679 1866.312
## 1996 2116.017 2137.537 2156.521 2173.546
## 1997 2297.142 2309.683 2322.567 2335.829
## 1998 2479.100 2488.579 2497.033 2505.433
## 1999 2630.912 2649.208 2669.550 2689.933
## 2000 NA NA NA NA# 绘制趋势拟合图
plot(x.fit$trend)
#查看随机波动(残差)值
x.fit$random## Jan Feb Mar Apr May Jun Jul
## 1993 NA NA NA NA NA NA 1.0088882
## 1994 0.9897286 0.9917699 1.0069996 1.0039379 1.0079296 1.0151591 0.9947887
## 1995 0.9982198 0.9571982 1.0006930 1.0123394 1.0122310 1.0092784 1.0128538
## 1996 0.9671360 1.0048115 1.0000891 1.0037414 1.0081543 1.0106040 0.9854906
## 1997 1.0006939 1.0088390 0.9998115 0.9993280 0.9918727 0.9988758 0.9963425
## 1998 1.0391783 0.9804972 0.9969239 0.9987216 0.9932308 0.9962666 1.0048103
## 1999 1.0145206 1.0128044 0.9892766 0.9843682 0.9776685 0.9823163 0.9878519
## 2000 0.9815488 1.0351058 0.9972322 0.9885893 0.9999390 0.9785257 NA
## Aug Sep Oct Nov Dec
## 1993 0.9880930 0.9748693 0.9584809 0.9492950 0.9863775
## 1994 1.0022841 1.0057935 0.9916461 1.0047607 1.0099323
## 1995 1.0064677 1.0009573 0.9948091 1.0044711 1.0085729
## 1996 0.9891157 1.0032647 0.9975563 1.0140385 1.0323984
## 1997 0.9935731 0.9933992 1.0090249 1.0092977 0.9718346
## 1998 1.0142266 1.0041262 1.0114724 1.0142280 0.9845558
## 1999 0.9972656 1.0086157 1.0280362 0.9949349 0.9973544
## 2000 NA NA NA NA NA# 绘制残差图
plot(x.fit$random)
1.2 Homework
1.2.1 1.
a5=scan('D:/New_Folder/Study_Programming/R_Programme/R-Time-Series/Homework/4th_homework/4.5.txt')
tsa5=ts(a5,start=1949)
ts.plot(tsa5)
该序列为显著的线性递增序列，故可以使用线性拟合或者holt两参数指数平滑法进行趋势拟合和预测。
t5=c(1949:2008)
a5.fit=lm(a5~t5)
summary(a5.fit)##
## Call:
## lm(formula = a5 ~ t5)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5320.0 -1142.8 730.5 1598.2 2328.6
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.771e+06 3.137e+04 -88.34 <2e-16 ***
## t5 1.449e+03 1.585e+01 91.38 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2127 on 58 degrees of freedom
## Multiple R-squared: 0.9931, Adjusted R-squared: 0.993
## F-statistic: 8351 on 1 and 58 DF, p-value: < 2.2e-16ts.plot(tsa5);abline(a5.fit,col=2)
故线性模型为 xt=−2.771∗106+1449t+ϵt,ϵt∼WN(0,21272)
a5.holt=HoltWinters(tsa5,gamma=F)
a5.holt## Holt-Winters exponential smoothing with trend and without seasonal component.
##
## Call:
## HoltWinters(x = tsa5, gamma = F)
##
## Smoothing parameters:
## alpha: 1
## beta : 1
## gamma: FALSE
##
## Coefficients:
## [,1]
## a 132802
## b 673plot(a5.holt)
# Forecast 10 years of population
library(forecast)
a5.fore=forecast(a5.holt,h=10)
plot(a5.fore)
1.2.2 2.
a6=scan('D:/New_Folder/Study_Programming/R_Programme/R-Time-Series/Homework/4th_homework/4.6.txt')
tsa6=ts(a6,start=1948,frequency=4)
ts.plot(tsa6)
该序列为显著的非线性递增序列，可以拟合二次型曲线、指数型曲线或其他曲线，也能使用holt两参数指数平滑法进行趋势拟合和预测。
考虑到为季度数据，我们不采用二次型曲线、指数型曲线或其他自变量单纯为t的曲线拟合。
下面采用holt两参数指数平滑法进行趋势拟合和预测。
a6.holt=HoltWinters(tsa6,gamma=F)
a6.holt## Holt-Winters exponential smoothing with trend and without seasonal component.
##
## Call:
## HoltWinters(x = tsa6, gamma = F)
##
## Smoothing parameters:
## alpha: 0.904549
## beta : 0.3212094
## gamma: FALSE
##
## Coefficients:
## [,1]
## a 5965.760
## b 146.505plot(a6.holt)
# Forecast 3 years of revenue
a6.fore=forecast(a6.holt,h=3)
plot(a6.fore)
1.2.3 3.
a7=scan('D:/New_Folder/Study_Programming/R_Programme/R-Time-Series/Homework/4th_homework/4.7.txt')
tsa7=ts(a7,frequency=12,start=c(1962,1))
ts.plot(tsa7)
观察时序图可知，该序列有显著趋势和周期效应。
下面尝试使用加法模型拟合该序列。
a7.fit1=decompose(tsa7,type='additive')
a7.fit1$figure## [1] -17.62196 -54.01259 33.38325 50.22700 110.14366 82.94054 28.62804
## [8] -14.84592 -53.10113 -49.72092 -75.54905 -40.47092plot(a7.fit1)
下面尝试用乘法模型拟合该序列
a7.fit2=decompose(tsa7,type='mult')
a7.fit2$figure## [1] 0.9752900 0.9239188 1.0481092 1.0718562 1.1567964 1.1178806 1.0416054
## [8] 0.9783800 0.9233615 0.9286178 0.8919580 0.9422262plot(a7.fit2)
选择乘法模型继续进行分析。
提取模型中的趋势项，并使用线性模型进行拟合。
ts.plot(a7.fit2$trend)
abline(lm(a7.fit2$trend~time(tsa7)),col="red")
reg<-lm(a7.fit2$trend~time(tsa7))
summary(reg)##
## Call:
## lm(formula = a7.fit2$trend ~ time(tsa7))
##
## Residuals:
## Min 1Q Median 3Q Max
## -7.802 -3.873 -1.056 3.644 10.955
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -4.269e+04 4.678e+02 -91.25 <2e-16 ***
## time(tsa7) 2.207e+01 2.379e-01 92.75 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.383 on 94 degrees of freedom
## (12 observations deleted due to missingness)
## Multiple R-squared: 0.9892, Adjusted R-squared: 0.9891
## F-statistic: 8603 on 1 and 94 DF, p-value: < 2.2e-16par(mfrow=c(2,2))
plot(reg,which=c(1:4))
显然回归模型拟合效果一般，残差出现了明显的自相关性。
使用该方法对模型进行预测需要提取出序列的季节因子，并与线性拟合的趋势项预测数据相乘。
首先我们提取出季节指数，得到季节指数sea，注意要先把results$seasonal项化为向量型，以便后面的计算。
如果前面分解采用加性模型，应该采用原始数据减去季节因素。
sea<-as.vector(a7.fit2$seasonal)
sea<-sea[1:12]
par(mfrow=c(1,1))
plot(sea,type = "l")
pre<-1971+(0:11)/12
pre<--42690.59+22.07*pre
pred<-pre*sea
cat('pred:',pred)## pred: 789.3802 749.5006 852.1739 873.453 944.798 915.07 854.5487 804.477 760.9361 766.9756 738.3376 781.6811最后得到pred项即为预测的结果。
1.2.4 4.
a8=scan('D:/New_Folder/Study_Programming/R_Programme/R-Time-Series/Homework/4th_homework/4.8.txt')
tsa8=ts(a8,frequency=12,start=c(1980,1))
ts.plot(tsa8)
此序列有曲线趋势,但没有固定周期效应,可在快速预测程序中用曲线拟合或曲线指数平滑进行预测。
此处用带季节的Holt-Winters三参数指数模型进行拟合及预测。
a8.fit=HoltWinters(tsa8)
a8.fit## Holt-Winters exponential smoothing with trend and additive seasonal component.
##
## Call:
## HoltWinters(x = tsa8)
##
## 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.48340plot(a8.fit)
a8.fore=forecast(a8.fit,h=3)
plot(a8.fore)