Sys.setlocale('LC_CTYPE', 'Chinese')
## [1] "Chinese (Simplified)_China.936"
load("D:\\Study_Programming\\R-documents\\Applied Statistics -codes and datas\\datas - Copy\\example\\ch11\\example11_1.RData")
example11_1=ts(example11_1,start = 2000)
par(mfrow=c(2,2),mai=c(0.6,0.6,0.12,0.1),cex=0.7)
plot(example11_1[,2],type = 'o',xlab='a.粮食产量序列',ylab = '粮食产量')
plot(example11_1[,3],type = 'o',xlab='b.居民消费水平序列',ylab = '居民消费水平')
plot(example11_1[,4],type = 'o',xlab='c.原煤产量序列',ylab = '原煤产量')
plot(example11_1[,5],type = 'o',xlab='d.CPI序列',ylab = 'CPI')
#指数平滑预测
#简单指数平滑预测
cpi_forecast=HoltWinters(example11_1[,5],beta = FALSE,gamma = FALSE)
cpi_forecast
## Holt-Winters exponential smoothing without trend and without seasonal component.
##
## Call:
## HoltWinters(x = example11_1[, 5], beta = FALSE, gamma = FALSE)
##
## Smoothing parameters:
## alpha: 0.2635414
## beta : FALSE
## gamma: FALSE
##
## Coefficients:
## [,1]
## a 102.3659
cpi_forecast$fitted#历史数据的拟合值
## Time Series:
## Start = 2001
## End = 2015
## Frequency = 1
## xhat level
## 2001 100.4000 100.4000
## 2002 100.4791 100.4791
## 2003 100.1420 100.1420
## 2004 100.4208 100.4208
## 2005 101.3377 101.3377
## 2006 101.4596 101.4596
## 2007 101.4702 101.4702
## 2008 102.3477 102.3477
## 2009 103.2839 103.2839
## 2010 102.2340 102.2340
## 2011 102.5149 102.5149
## 2012 103.2753 103.2753
## 2013 103.0973 103.0973
## 2014 102.9662 102.9662
## 2015 102.7116 102.7116
#绘制观测值、拟合值图
plot(example11_1[,5],type='o',xlab='时间',ylab='CPI')
lines(example11_1[,1][-1],cpi_forecast$fitted[,1],type = 'o',lty=2,col='blue')
legend(x='topleft',legend = c('观测值','拟合值'),lty=1:2,col=c(1,4))
#获取样本外的预测值
library(forecast)
cpi_forecast_1=forecast(cpi_forecast,h=1)
cpi_forecast_1
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 2016 102.3659 99.62907 105.1028 98.18025 106.5516
#实际值＆预测值
plot(cpi_forecast_1,type = 'o',xlab = '时间',ylab = 'CPI',main = '')
#Holt指数平滑预测
grain_forecast=HoltWinters(example11_1[,2],gamma = F)
grain_forecast
## Holt-Winters exponential smoothing with trend and without seasonal component.
##
## Call:
## HoltWinters(x = example11_1[, 2], gamma = F)
##
## Smoothing parameters:
## alpha: 0.666955
## beta : 0.5802066
## gamma: FALSE
##
## Coefficients:
## [,1]
## a 62211.953
## b 1108.538
#绘制拟合图
plot(example11_1[,2],type='o',xlab='时间',ylab='CPI')
lines(example11_1[,1][c(-1,-2)],grain_forecast$fitted[,1],type = 'o',lty=2,col='blue')
legend(x='topleft',legend = c('观测值','拟合值'),lty=1:2)
#预测
library(forecast)
grain_forecast_1=forecast(grain_forecast,h=1)
grain_forecast_1
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 2016 63320.49 61194.98 65446 60069.8 66571.18
#实际值、预测值图
plot(grain_forecast_1,type = 'o',xlab = '时间',ylab = '粮食产量',main = '')
#Winter指数平滑预测
load("D:\\Study_Programming\\R-documents\\Applied Statistics -codes and datas\\datas - Copy\\example\\ch11\\example11_4.RData")
example11_4=ts(example11_4[,3],start = 2010,frequency = 4)
example11_4
## Qtr1 Qtr2 Qtr3 Qtr4
## 2010 123 132 137 126
## 2011 130 138 142 132
## 2012 138 141 150 137
## 2013 143 147 158 143
## 2014 147 153 166 151
## 2015 159 163 174 161
#绘制各年销售量走势图、按年折叠的季度图、各年同季度图
library(forecast)
par(mfrow=c(1,3),mai=c(0.8,0.6,0.1,0.1),cex=0.7)
plot(example11_4,type='o',xlab='年份',ylab='销售量',pch=19,col=4,sub='a.各年销售量的走势图')
abline(v=c(2010:2015),lty=6,col='gray70')
seasonplot(example11_4,col=1,year.labels=T,xlab='季度',ylab='销售量',cex=0.8,main='',sub='b.销售量按年折叠的季度图')
monthplot(example11_4,xlab='季度',ylab='销售量',lty.base=2,col.base=2,sub='c.销售量的同季度图')
sale_forecast=HoltWinters(example11_4)
sale_forecast
## Holt-Winters exponential smoothing with trend and additive seasonal component.
##
## Call:
## HoltWinters(x = example11_4)
##
## Smoothing parameters:
## alpha: 0.3900487
## beta : 0.09609663
## gamma: 1
##
## Coefficients:
## [,1]
## a 165.6380871
## b 1.8531103
## s1 -0.9005488
## s2 1.2526248
## s3 10.8458341
## s4 -4.6380871
#Winter模型拟合图
plot(example11_4,type='o',xlab='时间',ylab='销售量')
lines(sale_forecast$fitted[,1],type = 'o',lty=2,col='blue')
legend(x='topleft',legend = c('观测值','拟合值'),lty=1:2,col=c(1:4))
#预测
library(forecast)
sale_forecast_1=forecast(sale_forecast,h=8)
sale_forecast_1
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 2016 Q1 166.5906 163.9802 169.2011 162.5984 170.5829
## 2016 Q2 170.5969 167.7580 173.4359 166.2551 174.9388
## 2016 Q3 182.0433 178.9557 185.1309 177.3212 186.7653
## 2016 Q4 168.4124 165.0578 171.7671 163.2819 173.5430
## 2017 Q1 174.0031 169.5015 178.5047 167.1185 180.8877
## 2017 Q2 178.0094 173.2621 182.7567 170.7490 185.2697
## 2017 Q3 189.4557 184.4443 194.4671 181.7915 197.1199
## 2017 Q4 175.8249 170.5320 181.1177 167.7302 183.9196
#预测图
plot(sale_forecast_1,type = 'o',xlab = '时间',ylab='销售量',main = '')
abline(v=2016,lty=6,col='gray')
#一元线性回归预测
fit=lm(粮食产量~年份,data=example11_1)
summary(fit)
##
## Call:
## lm(formula = 粮食产量 ~ 年份, data = example11_1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3465.1 -743.3 170.3 642.6 3463.6
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.478e+06 1.686e+05 -14.70 6.66e-10 ***
## 年份 1.260e+03 8.397e+01 15.01 5.05e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1548 on 14 degrees of freedom
## Multiple R-squared: 0.9415, Adjusted R-squared: 0.9373
## F-statistic: 225.2 on 1 and 14 DF, p-value: 5.053e-10
#各年预测值
predata=predict(fit,data.frame(年份=2000:2016))
predata
## 1 2 3 4 5 6 7 8
## 42753.92 44014.15 45274.38 46534.60 47794.83 49055.06 50315.29 51575.52
## 9 10 11 12 13 14 15 16
## 52835.75 54095.97 55356.20 56616.43 57876.66 59136.89 60397.12 61657.34
## 17
## 62917.57
#各年预测残差
res=residuals(fit)
res
## 1 2 3 4 5 6
## 3463.58162 1249.55324 431.42485 -3465.10353 -847.83191 -652.86029
## 7 8 9 10 11 12
## -511.08868 -1415.21706 35.15456 -1013.87382 -708.50221 504.46941
## 13 14 15 16
## 1081.34103 1056.91265 305.48426 486.55588
#观测值、预测值图
plot(2000:2016,predata,type = 'o',lty=2,col='blue',xlab='时间',ylab='粮食产量')
lines(2000:2015,example11_1[,2],type='o',pch=19)#按照课本代码example[,1]画不出来，需要修改，猜测是这一列非数值，未具体实验
legend(x='topleft',legend = c('观测值','预测值'),lty=1:2,col=c(1,4))
abline(v=2015,lty=6,col='gray')
#比较一元线性回归预测、Holt模型预测残差
residual_1=grain_forecast_1$residuals
residual_2=fit$residuals
plot(as.numeric(residual_1),ylim = c(-3000,3000),xlab = '',ylab = '残差',pch=1)
points(residual_2,pch=2,col='red')
abline(h=0)
legend(x='topright',legend = c('Holt模型预测残差','一元线性回归预测残差'),pch=1:2,col=c('black','red'),cex=0.8)
#分解预测
sale_compose=decompose(example11_4,type = 'multiplicative')
# names(sale_compose)
sale_compose$seasonal
## Qtr1 Qtr2 Qtr3 Qtr4
## 2010 0.9828110 1.0052503 1.0562235 0.9557153
## 2011 0.9828110 1.0052503 1.0562235 0.9557153
## 2012 0.9828110 1.0052503 1.0562235 0.9557153
## 2013 0.9828110 1.0052503 1.0562235 0.9557153
## 2014 0.9828110 1.0052503 1.0562235 0.9557153
## 2015 0.9828110 1.0052503 1.0562235 0.9557153
#绘制各成分图
plot(sale_compose,type='o')
#季节调整后的序列图
seasonal_adjust=example11_4/sale_compose$seasonal
plot(example11_4,xlab='时间',ylab='销售量',type='o',pch=19)
lines(seasonal_adjust,lty=2,type='o',col='blue')
legend(x='topleft',legend = c('销售量','销售量的季节调整'),lty=1:2)
#季节调整后的线性模型
x=1:24
fit_1=lm(seasonal_adjust~x)
fit_1
##
## Call:
## lm(formula = seasonal_adjust ~ x)
##
## Coefficients:
## (Intercept) x
## 124.312 1.692
#最终预测值
predata_1=predict(fit_1,data.frame(x=1:28))*rep(sale_compose$seasonal[1:4],7)
predata_1=ts(predata_1,start = 2010,frequency = 4)
predata_1
## Qtr1 Qtr2 Qtr3 Qtr4
## 2010 123.8384 128.3670 136.6635 125.2762
## 2011 130.4911 135.1716 143.8132 131.7455
## 2012 137.1438 141.9762 150.9628 138.2148
## 2013 143.7965 148.7808 158.1124 144.6841
## 2014 150.4492 155.5854 165.2620 151.1533
## 2015 157.1019 162.3899 172.4117 157.6226
## 2016 163.7546 169.1945 179.5613 164.0919
#预测残差
residual_3=example11_4-predict(fit_1,data.frame(x=1:24))*sale_compose$seasonal
residual_3=ts(residual_3,start = 2010,frequency = 4)
round(residual_3,4)
## Qtr1 Qtr2 Qtr3 Qtr4
## 2010 -0.8384 3.6330 0.3365 0.7238
## 2011 -0.4911 2.8284 -1.8132 0.2545
## 2012 0.8562 -0.9762 -0.9628 -1.2148
## 2013 -0.7965 -1.7808 -0.1124 -1.6841
## 2014 -3.4492 -2.5854 0.7380 -0.1533
## 2015 1.8981 0.6101 1.5883 3.3774
#实际值、预测值比较
plot(predata_1,type='o',lty=2,col='blue',xlab='时间',ylab='销售量')
lines(example11_4)
legend(x='topleft',legend=c('实际销售量','预测销售量'),lty=1:2,col=c("black" , "blue"),cex=0.7)
abline(v=2016,lty=6,col='grey')
#分解预测、Winter模型预测残差比较
sale_forecast_2=HoltWinters(example11_4)
residual_4=example11_4[-(1:4)]-sale_forecast_2$fitted[,1]
par(cex=0.7,mai=c(0.4,0.7,0.1,0.1))
plot(as.numeric(residual_3),xlab='',ylab = '残差',ylim = c(-5,5),pch=1)
abline(h=0)
points(as.numeric(residual_4),pch=2,col='red')
abline(h=0)
legend(x='topright',legend = c('分解预测的残差','Winter模型预测的残差'),pch=1:2,col=1:2,cex=0.8)
