原文链接:http://tecdat.cn/?p=3609
您要分析时间序列数据的第一件事就是将其读入R,并绘制时间序列。您可以使用scan()函数将数据读入R,该函数假定连续时间点的数据位于包含一列的简单文本文件中(点击文末“阅读原文”获取完整代码数据)。
读时间序列数据
数据集如下所示:
国王死亡年龄数据 60 43 67 50 56 42 50 65 68 43 65 34 ...
仅显示了文件的前几行。前三行包含对数据的一些注释,当我们将数据读入R时我们想要忽略它。我们可以通过使用scan()函数的“skip”参数来使用它,它指定了多少行。要忽略的文件顶部。要将文件读入R,忽略前三行,我们键入:
> kings [1] 60 43 67 50 56 42 50 65 68 43 65 34 47 34 49 41 13 35 53 56 16 43 69 59 48 [26] 59 86 55 68 51 33 49 67 77 81 67 71 81 68 70 77 56
在这种情况下,英国42位连续国王的死亡年龄已被读入变量“国王”。
一旦将时间序列数据读入R,下一步就是将数据存储在R中的时间序列对象中,这样就可以使用R的许多函数来分析时间序列数据。要将数据存储在时间序列对象中,我们使用R中的ts()函数。例如,要将数据存储在变量'kings'中作为R中的时间序列对象,我们键入:
Time Series: Start = 1 End = 42 Frequency = 1 [1] 60 43 67 50 56 42 50 65 68 43 65 34 47 34 49 41 13 35 53 56 16 43 69 59 48 [26] 59 86 55 68 51 33 49 67 77 81 67 71 81 68 70 77 56
您所拥有的时间序列数据集可能是以不到一年的固定间隔收集的,例如,每月或每季度。在这种情况下,您可以使用ts()函数中的'frequency'参数指定每年收集数据的次数。对于月度时间序列数据,您设置频率= 12,而对于季度时间序列数据,您设置频率= 4。
您还可以使用ts()函数中的“start”参数指定收集数据的第一年和该年度的第一个时间间隔。例如,如果第一个数据点对应于1986年第二季度,则设置start = c(1986,2)。
> birthstimeseries <- ts(births, frequency=12, start=c(1946,1))#转化成时间序列数据格式 > birthstimeseries Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec 1946 26.663 23.598 26.931 24.740 25.806 24.364 24.477 23.901 23.175 23.227 21.672 21.870 1947 21.439 21.089 23.709 21.669 21.752 20.761 23.479 23.824 23.105 23.110 21.759 22.073 1948 21.937 20.035 23.590 21.672 22.222 22.123 23.950 23.504 22.238 23.142 21.059 21.573 1949 21.548 20.000 22.424 20.615 21.761 22.874 24.104 23.748 23.262 22.907 21.519 22.025 1950 22.604 20.894 24.677 23.673 25.320 23.583 24.671 24.454 24.122 24.252 22.084 22.991 1951 23.287 23.049 25.076 24.037 24.430 24.667 26.451 25.618 25.014 25.110 22.964 23.981 1952 23.798 22.270 24.775 22.646 23.988 24.737 26.276 25.816 25.210 25.199 23.162 24.707 1953 24.364 22.644 25.565 24.062 25.431 24.635 27.009 26.606 26.268 26.462 25.246 25.180 1954 24.657 23.304 26.982 26.199 27.210 26.122 26.706 26.878 26.152 26.379 24.712 25.688 1955 24.990 24.239 26.721 23.475 24.767 26.219 28.361 28.599 27.914 27.784 25.693 26.881 1956 26.217 24.218 27.914 26.975 28.527 27.139 28.982 28.169 28.056 29.136 26.291 26.987 1957 26.589 24.848 27.543 26.896 28.878 27.390 28.065 28.141 29.048 28.484 26.634 27.735 1958 27.132 24.924 28.963 26.589 27.931 28.009 29.229 28.759 28.405 27.945 25.912 26.619 1959 26.076 25.286 27.660 25.951 26.398 25.565 28.865 30.000 29.261 29.012 26.992 27.897
同样, 1987年1月至1993年12月澳大利亚昆士兰州海滩度假小镇纪念品商店的月销售额(来自Wheelwright和Hyndman的原始数据, 1998)。我们可以通过输入以下内容将数据读入R:
Read 84 items > souvenirtimeseries <- ts(souvenir, frequency=12, start=c(1987,1)) > souvenirtimeseries Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec 1987 1664.81 2397.53 2840.71 3547.29 3752.96 3714.74 4349.61 3566.34 5021.82 6423.48 7600.60 19756.21 1988 2499.81 5198.24 7225.14 4806.03 5900.88 4951.34 6179.12 4752.15 5496.43 5835.10 12600.08 28541.72 1989 4717.02 5702.63 9957.58 5304.78 6492.43 6630.80 7349.62 8176.62 8573.17 9690.50 15151.84 34061.01 1990 5921.10 5814.58 12421.25 6369.77 7609.12 7224.75 8121.22 7979.25 8093.06 8476.70 17914.66 30114.41 1991 4826.64 6470.23 9638.77 8821.17 8722.37 10209.48 11276.55 12552.22 11637.39 13606.89 21822.11 45060.69 1992 7615.03 9849.69 14558.40 11587.33 9332.56 13082.09 16732.78 19888.61 23933.38 25391.35 36024.80 80721.71 1993 10243.24 11266.88 21826.84 17357.33 15997.79 18601.53 26155.15 28586.52 30505.41 30821.33 46634.38 104660.67
绘制时间序列
一旦你将时间序列读入R,下一步通常是制作时间序列数据的图,你可以用R中的plot.ts()函数做。
例如,为了绘制英国42位连续国王的死亡时间序列,我们输入:
> plot.ts(kingstimeseries)
我们可以从时间图中看出,可以使用加性模型来描述该时间序列,因为数据中的随机波动在大小上随时间大致恒定。
同样,为了绘制纽约市每月出生人数的时间序列,我们输入:
从这个时间序列我们可以看出,每月出生人数似乎有季节性变化:每年夏天都有一个高峰,每个冬天都有一个低谷。同样,似乎这个时间序列可能是用加性模型来描述的,因为季节性波动的大小随着时间的推移大致不变,似乎并不依赖于时间序列的水平,随机波动似乎也是随着时间的推移大小不变。
点击标题查阅往期内容
R语言中ARMA,ARIMA(Box-Jenkins),SARIMA和ARIMAX模型用于预测时间序列数据
01
02
03
04
同样,为了绘制澳大利亚昆士兰州海滩度假小镇纪念品商店每月销售的时间序列,我们输入:
在这种情况下,似乎加法模型不适合描述这个时间序列,因为季节性波动和随机波动的大小似乎随着时间序列的水平而增加。因此,我们可能需要转换时间序列以获得可以使用加法模型描述的变换时间序列。例如,我们可以通过计算原始数据的自然日志来转换时间序列:
> plot.ts(logsouvenirtimeseries)
在这里我们可以看到,对数变换时间序列中的季节性波动和随机波动的大小似乎随着时间的推移大致不变,并且不依赖于时间序列的水平。因此,可以使用加法模型来描述对数变换的时间序列。
分解时间序列
分解时间序列意味着将其分成其组成部分,这些组成部分通常是趋势分量和不规则分量,如果是季节性时间序列,则是季节性分量。
分解非季节性数据
非季节性时间序列由趋势分量和不规则分量组成。分解时间序列涉及尝试将时间序列分成这些分量,即估计趋势分量和不规则分量。
为了估计可以使用加性模型描述的非季节性时间序列的趋势分量,通常使用平滑方法,例如计算时间序列的简单移动平均值。
“TTR”R包中的SMA()函数可用于使用简单的移动平均值来平滑时间序列数据。要使用此功能,我们首先需要安装“TTR”R软件包 。一旦安装了“TTR”R软件包,就可以输入以下命令加载“TTR”R软件包:
然后,您可以使用“SMA()”功能来平滑时间序列数据。要使用SMA()函数,需要使用参数“n”指定简单移动平均值的顺序(跨度)。例如,要计算5阶的简单移动平均值,我们在SMA()函数中设置n = 5。
例如,如上所述,英国42位连续国王的死亡年龄的时间序列出现是非季节性的,并且可能使用加性模型来描述,因为数据中的随机波动大小基本上是恒定的。时间:
因此,我们可以尝试通过使用简单移动平均线进行平滑来估计此时间序列的趋势分量。要使用3阶简单移动平均值平滑时间序列,并绘制平滑时间序列数据,我们键入:
> kingstimeseriesSMA3 <- SMA(kingstimeseries,n=3) > plot.ts(kingstimeseriesSMA3)
在使用3阶简单移动平均值平滑的时间序列中,似乎存在相当多的随机波动。因此,为了更准确地估计趋势分量,我们可能希望尝试使用简单的移动平均值来平滑数据。更高阶。这需要一些试错,才能找到合适的平滑量。例如,我们可以尝试使用8阶的简单移动平均线:
> kingstimeseriesSMA8 <- SMA(kingstimeseries,n=8) > plot.ts(kingstimeseriesSMA8)
使用8阶简单移动平均值进行平滑的数据可以更清晰地显示趋势分量,我们可以看到英国国王的死亡年龄似乎已经从大约55岁降至大约38岁在最后的20位国王中,然后在第40位国王在时间序列的统治结束之后增加到大约73岁。
分解季节性数据
季节性时间序列由趋势组件,季节性组件和不规则组件组成。分解时间序列意味着将时间序列分成这三个组成部分:即估计这三个组成部分。
为了估计可以使用加性模型描述的季节性时间序列的趋势分量和季节性分量,我们可以使用R中的“decompose()”函数。该函数估计时间序列的趋势,季节和不规则分量。可以使用加性模型来描述。
函数“decompose()”返回一个列表对象作为结果,其中季节性组件,趋势组件和不规则组件的估计值存储在该列表对象的命名元素中,称为“季节性”,“趋势”和“随机” “ 分别。
例如,如上所述,纽约市每月出生人数的时间序列是季节性的,每年夏季和每年冬季都会出现高峰,并且可能使用加性模型来描述,因为季节性和随机波动似乎是随着时间的推移大小不变:
为了估计这个时间序列的趋势,季节性和不规则成分,我们输入:
> birthstimeseriescomponents <- decompose(birthstimeseries)
季节性,趋势和不规则成分的估计值现在存储在变量birthstimeseriescomponents seasonal,birthstimeseriescomponents seasonal,birthstimeseriescomponents trend和birthstimeseriescomponents $ random中。例如,我们可以通过键入以下内容输出季节性组件的估计值:
> birthstimeseriescomponents$seasonal # 得到季节成分 Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec 1946 -0.6771947 -2.0829607 0.8625232 -0.8016787 0.2516514 -0.1532556 1.4560457 1.1645938 0.6916162 0.7752444 -1.1097652 -0.3768197 1947 -0.6771947 -2.0829607 0.8625232 -0.8016787 0.2516514 -0.1532556 1.4560457 1.1645938 0.6916162 0.7752444 -1.1097652 -0.3768197 1948 -0.6771947 -2.0829607 0.8625232 -0.8016787 0.2516514 -0.1532556 1.4560457 1.1645938 0.6916162 0.7752444 -1.1097652 -0.3768197 1949 -0.6771947 -2.0829607 0.8625232 -0.8016787 0.2516514 -0.1532556 1.4560457 1.1645938 0.6916162 0.7752444 -1.1097652 -0.3768197 1950 -0.6771947 -2.0829607 0.8625232 -0.8016787 0.2516514 -0.1532556 1.4560457 1.1645938 0.6916162 0.7752444 -1.1097652 -0.3768197 1951 -0.6771947 -2.0829607 0.8625232 -0.8016787 0.2516514 -0.1532556 1.4560457 1.1645938 0.6916162 0.7752444 -1.1097652 -0.3768197 1952 -0.6771947 -2.0829607 0.8625232 -0.8016787 0.2516514 -0.1532556 1.4560457 1.1645938 0.6916162 0.7752444 -1.1097652 -0.3768197 1953 -0.6771947 -2.0829607 0.8625232 -0.8016787 0.2516514 -0.1532556 1.4560457 1.1645938 0.6916162 0.7752444 -1.1097652 -0.3768197 1954 -0.6771947 -2.0829607 0.8625232 -0.8016787 0.2516514 -0.1532556 1.4560457 1.1645938 0.6916162 0.7752444 -1.1097652 -0.3768197 1955 -0.6771947 -2.0829607 0.8625232 -0.8016787 0.2516514 -0.1532556 1.4560457 1.1645938 0.6916162 0.7752444 -1.1097652 -0.3768197 1956 -0.6771947 -2.0829607 0.8625232 -0.8016787 0.2516514 -0.1532556 1.4560457 1.1645938 0.6916162 0.7752444 -1.1097652 -0.3768197 1957 -0.6771947 -2.0829607 0.8625232 -0.8016787 0.2516514 -0.1532556 1.4560457 1.1645938 0.6916162 0.7752444 -1.1097652 -0.3768197 1958 -0.6771947 -2.0829607 0.8625232 -0.8016787 0.2516514 -0.1532556 1.4560457 1.1645938 0.6916162 0.7752444 -1.1097652 -0.3768197 1959 -0.6771947 -2.0829607 0.8625232 -0.8016787 0.2516514 -0.1532556 1.4560457 1.1645938 0.6916162 0.7752444 -1.1097652 -0.3768197
估计的季节性因素是在1月至12月期间给出的,并且每年都是相同的。最大的季节性因素是7月份(约1.46),最低的是2月份(约-2.08),表明7月出生率似乎达到高峰,2月出生低谷。
我们可以使用“plot()”函数绘制时间序列的估计趋势,季节和不规则分量,例如:
> plot(birthstimeseriescomponents)
上图显示了原始时间序列(顶部),估计趋势分量(从顶部开始的第二个),估计的季节性分量(从顶部开始的第三个)和估计的不规则分量(底部)。我们看到估计的趋势分量显示从1947年的大约24小幅下降到1948年的大约22小幅下降,随后从1959年开始稳步增加到大约27。
季节性调整
如果您有可以使用附加模型描述的季节性时间序列,则可以通过估计季节性成分来季节性地调整时间序列,并从原始时间序列中减去估计的季节性成分。我们可以使用“decompose()”函数计算的季节性成分的估计来做到这一点。
例如,要季节性调整纽约市每月出生人数的时间序列,我们可以使用“decompose()”估算季节性成分,然后从原始时间序列中减去季节性成分:
> birthstimeseriescomponents <- decompose(birthstimeseries) > birthstimeseriesseasonallyadjusted <- birthstimeseries - birthstimeseriescomponents$seasonal
然后我们可以使用“plot()”函数绘制经季节性调整的时间序列,输入:
> plot(birthstimeseriesseasonallyadjusted)
您可以看到季节性变化已从经季节性调整的时间序列中删除。经季节性调整的时间序列现在只包含趋势分量和不规则分量。
使用指数平滑的预测
指数平滑可用于对时间序列数据进行短期预测。
简单的指数平滑
如果您有一个时间序列可以使用具有恒定水平且没有季节性的附加模型来描述,则可以使用简单的指数平滑来进行短期预测。
简单指数平滑方法提供了一种估计当前时间点的水平的方法。平滑由参数alpha控制; 用于估计当前时间点的水平。alpha的值; α值接近于0意味着在对未来值进行预测时,最近的观察值很小。
Read 100 items > rainseries <- ts(rain,start=c(1813)) > plot.ts(rainseries)
你可以从图中看到大致恒定的水平(平均值保持恒定在25英寸左右)。随着时间的推移,时间序列中的随机波动似乎大致不变,因此使用加性模型描述数据可能是合适的。因此,我们可以使用简单的指数平滑进行预测。
为了使用R中的简单指数平滑进行预测,我们可以使用R中的“HoltWinters()”函数拟合一个简单的指数平滑预测模型。要使用HoltWinters()进行简单的指数平滑,我们需要设置参数beta = FALSE和HoltWinters()函数中的gamma = FALSE(β和gamma参数用于Holt的指数平滑,或Holt-Winters指数平滑,如下所述)。
HoltWinters()函数返回一个列表变量,该变量包含多个命名元素。
例如,要使用简单的指数平滑来预测伦敦年降雨量的时间序列,我们输入:
> rainseriesforecasts <- HoltWinters(rainseries, beta=FALSE, gamma=FALSE) > rainseriesforecasts Smoothing parameters: alpha: 0.02412151 beta : FALSE gamma: FALSE Coefficients: [,1] a 24.67819
HoltWinters()的输出告诉我们alpha参数的估计值约为0.024。这非常接近零,告诉我们预测是基于最近和最近的观察结果(虽然对最近的观察更加重视)。
默认情况下,HoltWinters()仅对我们原始时间序列所涵盖的相同时间段进行预测。在这种情况下,我们的原始时间序列包括1813年至1912年伦敦的降雨量,所以预测也是1813年至1912年。
在上面的例子中,我们将HoltWinters()函数的输出存储在列表变量“rainseriesforecasts”中。HoltWinters()的预测存储在这个名为“fits”的列表变量的命名元素中,因此我们可以通过输入以下内容来获取它们的值:
> rainseriesforecasts$fitted Time Series: Start = 1814 End = 1912 Frequency = 1 xhat level 1814 23.56000 23.56000 1815 23.62054 23.62054 1816 23.57808 23.57808 1817 23.76290 23.76290 1818 23.76017 23.76017 1819 23.76306 23.76306 1820 23.82691 23.82691 ... 1905 24.62852 24.62852 1906 24.58852 24.58852 1907 24.58059 24.58059 1908 24.54271 24.54271 1909 24.52166 24.52166 1910 24.57541 24.57541 1911 24.59433 24.59433 1912 24.59905 24.59905
我们可以通过键入以下内容来绘制原始时间序列与预测:
> plot(rainseriesforecasts)
该图显示原始时间序列为黑色,预测显示为红线。预测的时间序列比原始数据的时间序列要平滑得多。
作为预测准确性的度量,我们可以计算样本内预测误差的平方误差之和,即我们原始时间序列所涵盖的时间段的预测误差。平方误差之和存储在名为“SSE”的列表变量“rainseriesforecasts”的命名元素中,因此我们可以通过键入以下内容来获取其值:
> rainseriesforecasts$SSE [1] 1828.855
也就是说,这里的平方误差之和为1828.855。
在简单的指数平滑中,通常使用时间序列中的第一个值作为级别的初始值。例如,在伦敦的降雨时间序列中,1813年降雨量的第一个值为23.56(英寸)。您可以使用“l.start”参数指定HoltWinters()函数中水平的初始值。例如,要将级别的初始值设置为23.56进行预测,我们键入:
> HoltWinters(rainseries, beta=FALSE, gamma=FALSE, l.start=23.56)
如上所述,默认情况下,HoltWinters()仅对原始数据所涵盖的时间段进行预测,即降雨时间序列为1813-1912。我们可以使用R“forecast”包中的“forecast.HoltWinters()”函数对更多时间点进行预测。要使用forecast.HoltWinters()函数,我们首先需要安装“预测”R包(有关如何安装R包的说明,请参阅如何安装R包)。
安装“预测”R软件包后,您可以键入以下命令加载“预测”R软件包:
> library("forecast")
当使用forecast.HoltWinters()函数作为其第一个参数(输入)时,您将使用HoltWinters()函数传递给您已经拟合的预测模型。例如,在降雨时间序列的情况下,我们将使用HoltWinters()的预测模型存储在变量“rainseriesforecasts”中。您可以使用forecast.HoltWinters()中的“h”参数指定要进行预测的其他时间点数。例如,要使用forecast.HoltWinters()预测1814-1820(8年以上)的降雨量,我们输入:
> rainseriesforecasts2 <- forecast.HoltWinters(rainseriesforecasts, h=8) > rainseriesforecasts2 Point Forecast Lo 80 Hi 80 Lo 95 Hi 95 1913 24.67819 19.17493 30.18145 16.26169 33.09470 1914 24.67819 19.17333 30.18305 16.25924 33.09715 1915 24.67819 19.17173 30.18465 16.25679 33.09960 1916 24.67819 19.17013 30.18625 16.25434 33.10204 1917 24.67819 19.16853 30.18785 16.25190 33.10449 1918 24.67819 19.16694 30.18945 16.24945 33.10694 1919 24.67819 19.16534 30.19105 16.24701 33.10938 1920 24.67819 19.16374 30.19265 16.24456 33.11182
forecast.HoltWinters()函数为您提供一年的预测,预测的预测间隔为80%,预测的预测间隔为95%。例如,1920年的预测降雨量约为24.68英寸,95%的预测间隔为(16.24,33.11)。
要绘制forecast.HoltWinters()所做的预测,我们可以使用“plot.forecast()”函数:
> plot.forecast(rainseriesforecasts2)
使用R语言进行时间序列(arima,指数平滑)分析(下):https://developer.aliyun.com/article/1493894