[读者文库]Introduction to Time Series Analysis and Forecasting(using R)

Example 4.2
The US CPI data are in the second column of the array called cpi.data in which the first column is the month of the year. For this case we use the firstsmooth function twice to obtain the double exponential smoothing as

1. cpi.smooth1<-firstsmooth(y=cpi.data[,2],lambda=0.3)cpi.smooth2<-firstsmooth(y=cpi.smooth1,lambda=0.3)cpi.hat<-2*cpi.smooth1-cpi.smooth2 #Equation 4.23plot(cpi.data[,2],type="p", pch=16,cex=.5,xlab='Date',ylab='CPI',xaxt='n')axis(1, seq(1,120,24), cpi.data[seq(1,120,24),1])lines(cpi.hat)

复制代码

Note that the fitted values are obtained using Eq. (4.23). Also the corresponding command using Holt–Winters function is

1. HoltWinters(cpi.data[,2],alpha=0.3, beta=0.3, gamma=FALSE)

复制代码

Example 4.4

The average speed data are in the second column of the array called speed.data in which the first column is the index for the week.


To find the “best” smoothing constant, we will use the firstsmooth function for various lambda values and obtain the sum of squared one-step-ahead prediction error (SSE) for each. The lambda value that minimizes the sum of squared prediction errors is deemed the “best” lambda. The obvious option is to apply firstsmooth function in a for loop to obtain SSE for various lambda values. Even though in this case this may not be an issue, in many cases for loops can slow down the computations in R and are to be avoided if possible. We will do that using sapply function.

1. lambda.vec<-seq(0.1, 0.9, 0.1)
   
2. sse.speed<-function(sc){measacc.fs(speed.data[1:78,2],sc)[1]}
   
3. sse.vec<-sapply(lambda.vec, sse.speed)
   
4. opt.lambda<-lambda.vec[sse.vec == min(sse.vec)]
   
5. plot(lambda.vec, sse.vec, type="b", main = "SSE vs. lambda\n",
   
6. xlab='lambda\n',ylab='SSE')
   
7. abline(v=opt.lambda, col = 'red')
   
8. mtext(text = paste("SSE min = ", round(min(sse.vec),2), "\n lambda
   
9. = ", opt.lambda))

复制代码

Example 4.5
We will first try to find the best lambda for the CPI data using first-order exponential smoothing.We will also plot ACF of the data. Note that we will use the data up to December 2003.

1. lambda.vec<-c(seq(0.1, 0.9, 0.1), .95, .99)sse.cpi<-function(sc){measacc.fs(cpi.data[1:108,2],sc)[1]}sse.vec<-sapply(lambda.vec, sse.cpi)opt.lambda<-lambda.vec[sse.vec == min(sse.vec)]plot(lambda.vec, sse.vec, type="b", main = "SSE vs. lambda\n",xlab='lambda\n',ylab='SSE', pch=16,cex=.5)acf(cpi.data[1:108,2],lag.max=25

复制代码

Example 4.6
The function for the Trigg–Leach smoother is given as:

1. tlsmooth<-function(y,gamma,y.tilde.start=y[1],lambda.start=1){
   
2. T<-length(y)
   
3. #Initialize the vectors
   
4. Qt<-vector()
   
5. Dt<-vector()
   
6. y.tilde<-vector()
   
7. lambda<-vector()
   
8. err<-vector()
   
9. #Set the starting values for the vectors
   
10. lambda[1]=lambda.start
    
11. y.tilde[1]=y.tilde.start
    
12. Qt[1]<-0
    
13. Dt[1]<-0
    
14. err[1]<-0
    
15. for (i in 2:T){
    
16. err[i]<-y[i]-y.tilde[i-1]
    
17. Qt[i]<-gamma*err[i]+(1-gamma)*Qt[i-1]
    
18. Dt[i]<-gamma*abs(err[i])+(1-gamma)*Dt[i-1]
    
19. lambda[i]<-abs(Qt[i]/Dt[i])
    
20. y.tilde[i]=lambda[i]*y[i] + (1-lambda[i])*y.tilde[i-1]
    
21. }
    
22. return(cbind(y.tilde,lambda,err,Qt,Dt))
    
23. }
    
24. #Obtain the TL smoother for Dow Jones Index
    
25. out.tl.dji<-tlsmooth(dji.data[,2],0.3)

复制代码

Example 4.7
The clothing sales data are in the second column of the array called closales.data in which the first column is the month of the year. We will use the data up to December 2002 to fit the model and make forecasts for the coming year (2003). We will use Holt–Winters function given in stats package. The model is additive seasonal model with all parameters equal to 0.2.

1. dat.ts = ts(closales.data[,2], start = c(1992,1), freq = 12)
   
2. y1<-closales.data[1:132,]
   
3. # convert data to ts object
   
4. y1.ts<-ts(y1[,2], start = c(1992,1), freq = 12)
   
5. clo.hw1<-HoltWinters(y1.ts,alpha=0.2,beta=0.2,gamma=0.2,seasonal
   
6. ="additive")
   
7. plot(y1.ts,type="p", pch=16,cex=.5,xlab='Date',ylab='Sales')
   
8. lines(clo.hw1$fitted[,1])
   
9. #Forecast the the sales for 2003
   
10. y2<-closales.data[133:144,]
    
11. y2.ts<-ts(y2[,2],start=c(2003,1),freq=12)
    
12. y2.forecast<-predict(clo.hw1, n.ahead=12, prediction.interval
    
13. = TRUE)
    
14. plot(y1.ts,type="p", pch=16,cex=.5,xlab='Date',ylab='Sales',
    
15. xlim=c(1992,2004))
    
16. points(y2.ts)
    
17. lines(y2.forecast[,1])
    
18. lines(y2.forecast[,2])
    
19. lines(y2.forecast[,3])

复制代码