通过R语言实现平稳时间序列的建模 |
您所在的位置:网站首页 › r语言建立var模型 › 通过R语言实现平稳时间序列的建模 |
目录 1. 建模流程 2. 序列平稳性检验和纯随机性检验 2.1 图检验 2.2 单位根检验 3. 模型选择 4. 参数估计 5. 模型检验 5.1 模型显著性检验 5.2 参数显著性检验 6. 模型优化 6.1 AIC准则 6.2 BIC准则 7. 预测 1. 建模流程1.1 序列平稳性检验+纯随机性检验 1.2 模型选择 1.3 参数估计 1.4 模型检验 1.5 模型优化 1.6 预测 2. 序列平稳性检验和纯随机性检验 2.1 图检验主要适用于趋势或者周期比较明显的序列,具有一定主观性。 (1)绘制时序图与自相关图: x Box.test(dwelling,lag=12,type="Ljung-Box") #延迟12阶 Box-Ljung test data: dwelling X-squared = 329.2, df = 12, p-value < 2.2e-16白噪声检验显示延迟为6、12阶时p值均远远小于显著性水平0.05,因此显著拒绝序列为纯随机系列的原假设。因此判断该序列为非纯随机序列。 2.2 单位根检验2.2.1 DF检验 仅适用于最简单的1阶模型AR(1)的平稳性检验。 调用aTSA中的的adf.test 函数进行DF检验: install.packages("aTSA") library(aTSA) adf.test(dwelling,nlag=2) #nlag:最高延迟阶数,若nlag=1即输出自回归0阶延迟平稳性检验结果;若nlag=2即输出自回归0至1阶延迟平稳性检验结果。 #输出结果 > adf.test(dwelling,nlag=2) Augmented Dickey-Fuller Test alternative: stationary Type 1: no drift no trend lag ADF p.value [1,] 0 -2.46 0.0162 [2,] 1 -1.90 0.0573 Type 2: with drift no trend lag ADF p.value [1,] 0 -4.81 0.01 [2,] 1 -3.68 0.01 Type 3: with drift and trend lag ADF p.value [1,] 0 -4.97 0.0100 [2,] 1 -3.89 0.0185 ---- Note: in fact, p.value = 0.01 means p.value adf.test(dwelling) Augmented Dickey-Fuller Test alternative: stationary Type 1: no drift no trend lag ADF p.value [1,] 0 -2.46 0.0162 [2,] 1 -1.90 0.0573 [3,] 2 -1.38 0.1832 [4,] 3 -1.10 0.2827 Type 2: with drift no trend lag ADF p.value [1,] 0 -4.81 0.010 [2,] 1 -3.68 0.010 [3,] 2 -2.43 0.161 [4,] 3 -1.94 0.350 Type 3: with drift and trend lag ADF p.value [1,] 0 -4.97 0.0100 [2,] 1 -3.89 0.0185 [3,] 2 -2.60 0.3254 [4,] 3 -2.03 0.5579 ---- Note: in fact, p.value = 0.01 means p.value arima(dwelling,order=c(1,0,0),method="ML") Call: arima(x = dwelling, order = c(1, 0, 0), method = "ML") Coefficients: ar1 intercept 0.5929 40.5062 s.e. 0.0877 4.6893 sigma^2 estimated as 331.3: log likelihood = -380.41, aic = 766.82第一行输出的是参数估计值,第二行输出的是参数估计值的样本标准差。 根据输出结果,得到模型: 或者,通过输出的参数得到 因此得到AR(1)模型为: 如果残差序列为非白噪声检验,则意味着残差序列中依旧含有未被模型提取的相关信息,这就说明这个拟合模型不够有效。因此,模型的显著性检验即残差序列的白噪声检验。 原假设为该序列的残差均为0。 fit=arima(dwelling,order=c(1,0,0),method="ML") #fit为arima函数拟合结果 ts.diag(fit) #进行模型显著性检验,输出四个图 #左上的图为残差序列自相关图;右上的图为残差序列偏相关图 #左下的图为残差序列白噪声检验图;右下为残差序列正态性检验QQ图输出四图如下: 左下方的残差序列自相关图横轴为延迟阶数,纵轴为该延阶数纯随机性检验(Q统计量)的p值。若所有的点都在虚线上方,即所有Q统计量的p值均大于显著性水平0.05,此时可认为该模型显著成立。 右下方的残差序列正态性检验QQ图横轴为正态分布的分位数,纵轴为样本分位数,若图上的点密集分布在对称线左右,则可以认为该序列服从正态分布。该图可用以对序列进行正态分布假定的检验。 如图,由于残差序列自相关图上所有点均在虚线上方,即所有p值均显著大于0.05,可以认为该模型的残差序列为白噪声序列,即该模型通过了显著性检验。 5.2 参数显著性检验即检验未知参数是否显著非零。若某个参数不显著非零,则可以从拟合模型中剔除,使得模型精简。 5.2.1 法一 当样本较大时,序列近似于正态分布,因此只要参数估计值大于其两倍标准差,我们就可以认为该参数显著非零。 此时只需要再使用一次arima函数: arima(dwelling,order=c(1,0,0),method="ML") #输出结果 > arima(dwelling,order=c(1,0,0),method="ML") Call: arima(x = dwelling, order = c(1, 0, 0), method = "ML") Coefficients: ar1 intercept 0.5929 40.5062 s.e. 0.0877 4.6893因为0.5929>0.0877*2,20.5062>4.6893*2,因此两个参数均显著非零。 5.2.2 法二 先构造参数显著性检验的t统计量,再通过pt函数求得p值: t=abs(fit$coef)/sqrt(diag(fit$var.coef)) #构造t统计量 pt(t,length(dwelling)-length(fit$coef),lower.tail = F) #输出p值 #输出结果 > t=abs(fit$coef)/sqrt(diag(fit$var.coef)) > pt(t,length(dwelling)-length(fit$coef),lower.tail = F) ar1 intercept 7.828497e-10 1.343084e-13两个参数t统计量的p值均远小于显著性水平0.05,因此两个参数均显著非零。 6. 模型优化一个好的模型是拟合精度与未知参数个数的综合最优配置,我们可以通过最小信息量检查AIC以及它的改进BIC准则来选出相对的最优模型。 6.1 AIC准则 AIC(fit) #输出结果 AIC(fit) [1] 766.8243可以对比多个模型的AIC值,使得AIC达到最小值的模型为相对最优模型。 6.2 BIC准则 BIC(fit) #输出结果 > BIC(fit) [1] 774.2563可以对比多个模型的BIC值,使得BIC达到最小值的模型为相对最优模型。 也可以同时考虑AIC函数值与BIC函数值,使得二者达到最小的模型为相对最优模型。 6.3 auto.arima 函数 该函数基于信息量最小原则,可以在一定范围内自动识别最优模型的阶数,可以起到简化模型优化问题的作用。 install.packages("forecast") library("forecast") auto.arima(dwelling,max.p=5,max.q=5,ic="bic") #max.p:指定自相关系数的最高阶数,默认为5; #max.q:指定移动平均系数的最高阶数,默认为5; #ic:指定信息量准则,默认为ic="aic",即AIC准则;ic="bic"即BIC准则;此外还有ic="aicc"。 #输出结果 > auto.arima(dwelling,max.p=5,max.q=5,ic="bic") Series: dwelling ARIMA(1,0,1)(1,0,0)[12] with zero mean Coefficients: ar1 ma1 sar1 0.9648 -0.4870 0.3971 s.e. 0.0297 0.0951 0.1375 sigma^2 estimated as 296.6: log likelihood=-375.97 AIC=759.95 AICc=760.43 BIC=769.86考虑到输出的第二个参数值-0.4870 |
今日新闻 |
推荐新闻 |
CopyRight 2018-2019 办公设备维修网 版权所有 豫ICP备15022753号-3 |