R语言时间序列代码整理

最新推荐文章于 2023-01-13 20:09:15 发布

原创最新推荐文章于 2023-01-13 20:09:15 发布 · 5.7k 阅读

50 ·

CC 4.0 BY-SA版权

数据挖掘专栏收录该内容

1 篇文章

订阅专栏

本文整理了R语言中时间序列分析的关键步骤，包括读取数据、平滑化、季节性分解、指数平滑预测（ETS）和ARIMA模型。通过示例代码展示了如何进行简单移动平均、STL分解、ETS模型拟合以及ARIMA模型的选择和预测，为时间序列预测提供了实用的R代码参考。

R语言时间序列代码整理

时间序列：
data<-read.csv('/Users/ji_fanyang/Desktop/test2.csv')
myseries<- ts(data[1],start=c(2017,1),end=c(2020,7),frequency=12)
install.packages('forecast')
library(forecast)

1.时序的平滑化
平滑曲线，居中移动平均
opar<-par(no.readonly=TRUE)
par(mfrow=c(2,2))
ylim<-c(min(myseries),max(myseries))
plot(myseries,main="Raw time series")
plot(ma(myseries,3),main="simple moving average(k=3)",ylim=ylim)
plot(ma(myseries,7),main="simple moving average(k=7)",ylim=ylim)
plot(ma(myseries,15),main="simple moving average(k=15)",ylim=ylim)

2.季节性分解（趋势因子，季节因子，随机因子）
stl=(ts,s.window=,t.window=)，为相加模型

S.window:控制季节效应变化的速度，s.window=“periodic”,使得季节效应在每个月都是一致的
t.window:控制趋势变化的速度

myseries2<-log(myseries)
fit<-stl(myseries,s.window="period")
plot(fit)
fit$time.series
exp(fit$time.series):将结果转化为原始尺寸

将季节分解进行可视化
par(mfrow=c(2,1))
monthplot(myseries,xlab='',ylab='')
seasonplot(myseries,year.labels="TRUE")

3.指数预测模型
fit<-ets(log(myseries),model='AAA')
fit(输出）
ETS(A,A,A)

Call:
ets(y = log(myseries), model = "AAA")

Smoothing parameters:
    alpha = 1e-04 #水平项
    beta = 1e-04 #趋势项，小意味着近期观测值的斜率不需要更新
    gamma = 4e-04 #季节项

Initial states:
    l = 6.2117
    b = -0.0077
    s = -0.5525 -0.056 -0.7451 -1.2488 1.635 0.4444
           1.7078 0.8398 -1.1126 0.2562 -0.2587 -0.9094

sigma: 3.3361

AIC AICc BIC
279.3346 303.8146 309.2750

accuracy(fit) #
ME RMSE MAE MPE MAPE MASE
Training set 0.03006177 2.643563 2.103576 -33.22007 56.21131 0.6137234

pred <-forecast(fit,5) #预测接下来5个月

plot(pred,main='',ylab='log(myseries)',xlab='Time')
pred$mean <-exp(pred$mean)    #返回均值
pred$lower <-exp(pred$lower)
pred$upper <-exp(pred$upper)   #分别包含80%和95%置信区间的下界和上界
p<-cbind(pred$mean,pred$lower,pred$upper)

3.ets()和自动预测
fit2<-ets(myseries)
生成数据：
fit2
ETS(M,N,M)

Call:
ets(y = myseries)

Smoothing parameters:
alpha = 0.1207
gamma = 0.1438

Initial states:
    l = 88818.7326
    s = 0.6389 0.0308 1.2384 0.0096 1.3914 2.5175
           0.4078 3.1015 0.013 0.6068 1.7915 0.2527

sigma: 1.1362

AIC AICc BIC
971.3118 989.0896 997.7298
> pred2<-forecast(fit2)

4.ARIMA模型
1.验证序列的平稳性
library(forecast)
library(tseries)
ndiff(myseries) #找到最优值
dmyseries <-diff(myseries,1)
adf.test(dmyseries)

2.选择模型
Acf(dmyseries)
Pacf(dmyseries)

3.拟合模型
fit<-arima(dmyseries,order=c(0,1,1))

4.模型评价（QQ图）
qqnorm(fit$residuals)
qqline(fit$residuals) #如果数据满足正态分布，则数据中的点会落在图中的线上
Box.test(fit$residuals,type=“Ljung-Box”) #可以检验残差的自相关系数是否都为0

5.预测：如果残差满足QQ图或者自相关系数为0
forecast(fit,3)
accuracy(forecast)

—Arima自动预测
fit<-auto.arima(myseries)