语言时间序列年月日_R语言系列 时间序列分析

【免责声明：本文用于教学】

时间序列分析

基础操作

数据输入

d <- c(10,15,10,10,12,10,7,7,10,14,8,17,14,18,3,9,11,10,6,12,14,10,25,29,33,33,12,19,16,19,19,12,34,15,36,29,26,21,17,19,13,20,24,12,6,14,6,12,9,11,17,12,8,14,14,12,5,8,10,3,16,8,8,7,12,6,10,8,10,5)
##convert it to a time series
is.ts(d)

## [1] FALSE

d <- ts(d, start=min(201401))
d

## Time Series:
## Start = 201401
## End = 201470
## Frequency = 1
## [1] 10 15 10 10 12 10 7 7 10 14 8 17 14 18 3 9 11 10 6 12 14 10 25 29 33
## [26] 33 12 19 16 19 19 12 34 15 36 29 26 21 17 19 13 20 24 12 6 14 6 12 9 11
## [51] 17 12 8 14 14 12 5 8 10 3 16 8 8 7 12 6 10 8 10 5

plot.ts(d,ylab="y", xlab="day")

ACF和PACF

par(mfrow=c(1,2))
##autocovariance function ACF 自相关
acf(d)
acf(d,plot=F)

##
## Autocorrelations of series 'd', by lag
##
## 0 1 2 3 4 5 6 7 8 9 10
## 1.000 0.506 0.539 0.374 0.291 0.258 0.148 0.270 0.186 0.178 0.258
## 11 12 13 14 15 16 17 18
## 0.207 0.226 0.138 -0.027 -0.053 -0.112 -0.139 -0.155

##partial autocovariance function PACF 偏自相关
pacf(d)

pacf(d,plot=F)

##
## Partial autocorrelations of series 'd', by lag
##
## 1 2 3 4 5 6 7 8 9 10 11
## 0.506 0.380 0.020 -0.055 0.046 -0.062 0.207 0.032 -0.069 0.163 0.037
## 12 13 14 15 16 17 18
## -0.023 -0.055 -0.318 -0.101 0.087 -0.091 -0.060

检验

单位根检验 : 非平稳 vs : 平稳

library(fUnitRoots)

## Loading required package: timeDate

## Loading required package: timeSeries

## Loading required package: fBasics

unitrootTest(d)

##
## Title:
## Augmented Dickey-Fuller Test
##
## Test Results:
## PARAMETER:
## Lag Order: 1
## STATISTIC:
## DF: -1.2139
## P VALUE:
## t: 0.2041
## n: 0.4322
##
## Description:
## Mon Dec 7 00:41:14 2020 by user:

纯随机性 :白噪声 vs : 非白噪声

Box.test(d,type="Ljung-Box")

##
## Box-Ljung test
##
## data: d
## X-squared = 18.7, df = 1, p-value = 1.53e-05

Box.test(d,type="Ljung-Box",lag=6)

##
## Box-Ljung test
##
## data: d
## X-squared = 64.016, df = 6, p-value = 6.85e-12

Box.test(d,type="Ljung-Box",lag=12)

##
## Box-Ljung test
##
## data: d
## X-squared = 88.975, df = 12, p-value = 7.794e-14

做差分

d1 <- diff(d)
plot.ts(d1,ylab="diff", xlab="day")

par(mfrow=c(1,2))
##autocovariance function ACF 自相关
acf(d1)
acf(d1,plot=F)

##
## Autocorrelations of series 'd1', by lag
##
## 0 1 2 3 4 5 6 7 8 9 10
## 1.000 -0.529 0.195 -0.080 -0.059 0.092 -0.256 0.216 -0.075 -0.070 0.101
## 11 12 13 14 15 16 17 18
## -0.048 0.104 0.075 -0.142 0.045 -0.032 -0.026 -0.022

##partial autocovariance function PACF 偏自相关
pacf(d1)

pacf(d1,plot=F)

##
## Partial autocorrelations of series 'd1', by lag
##
## 1 2 3 4 5 6 7 8 9 10 11
## -0.529 -0.119 -0.039 -0.142 -0.011 -0.283 -0.095 0.015 -0.175 -0.091 -0.015
## 12 13 14 15 16 17 18
## 0.024 0.272 0.038 -0.121 0.043 0.002 -0.057

##单位根检验 H0: 非平稳
library(fUnitRoots)
unitrootTest(d1)

##
## Title:
## Augmented Dickey-Fuller Test
##
## Test Results:
## PARAMETER:
## Lag Order: 1
## STATISTIC:
## DF: -7.9309
## P VALUE:
## t: < 2.2e-16
## n: 0.0461
##
## Description:
## Mon Dec 7 00:41:14 2020 by user:

##纯随机性 H0:白噪声
Box.test(d1,type="Ljung-Box")

##
## Box-Ljung test
##
## data: d1
## X-squared = 20.195, df = 1, p-value = 6.993e-06

Box.test(d1,type="Ljung-Box",lag=6)

##
## Box-Ljung test
##
## data: d1
## X-squared = 29.458, df = 6, p-value = 4.982e-05

Box.test(d1,type="Ljung-Box",lag=12)

##
## Box-Ljung test
##
## data: d1
## X-squared = 35.943, df = 12, p-value = 0.0003308

ts函数

R中自带的ts()函数可以直接生成时序数据，调用公式如下：

ts(data = NA, start = , end = , frequency = , ...)
frequency：取值包括1,4,12,分别代表年度,季度,月度数据

##生成季度数据,且以1959年第二季度为起始点(月度数据同理)
ts(1:10, frequency = 4, start = c(1959, 2))

## Qtr1 Qtr2 Qtr3 Qtr4
## 1959 1 2 3
## 1960 4 5 6 7
## 1961 8 9 10

# 同时生成多条时间序列
data <- ts(matrix(rnorm(180), 60, 3), start = c(2012, 1), frequency = 12,names = c('a','b','c'))
head(data)

## a b c
## [1,] 1.38880495 -1.7780026 -0.8546904
## [2,] -0.03141699 0.9623493 1.2816379
## [3,] 0.69086949 0.5687210 -2.0918558
## [4,] 0.86851666 -0.9937985 2.2763149
## [5,] -0.58213129 -0.2361856 0.2417787
## [6,] 0.69791200 -1.3831174 0.6326396

plot(data)

ts.data <- ts(data,start=c(1996,10),freq=12)
plot(ts.data)

# 使用window函数选取部分数据
data_subset <- window(ts.data,start=c(1998,1),end=c(1998,6))
data_subset

## a b c
## Jan 1998 0.09813527 0.2172911 0.10958001
## Feb 1998 1.78146726 -0.7813561 0.39675998
## Mar 1998 -1.07601484 0.2697066 -0.87731481
## Apr 1998 0.53832758 1.5250992 -0.09629001
## May 1998 0.71610457 0.8583480 -0.29970496
## Jun 1998 -0.33019984 0.2186358 -1.36307364

时间序列分解

简单移动平均

R中的forecast包的ma()函数、zoo包的rollmean()函数和TTR包的SMA()函数均可实现简单移动平均处理，这里我们选用TTR包的SMA函数

# nhtemp数据集为康涅狄格州纽黑文地区从1912年至1971年每年的平均气温时序数据
data(nhtemp)
ylim <- c(45,55)
library(TTR)

##
## Attaching package: 'TTR'

## The following object is masked from 'package:fBasics':
##
## volatility

par(mfrow=c(2,2))
plot(nhtemp, main="Raw time series")
plot(SMA(nhtemp,3), main="SMA(k=3)",ylim=ylim)
plot(SMA(nhtemp,5), main="SMA(k=5)",ylim=ylim)
plot(SMA(nhtemp,7), main="SMA(k=7)",ylim=ylim)

stl函数

stl(ts, s.window=, t.window=)

s.window：控制季节效应变化的速度,设定为"periodic"时可使得季节效应在各年间都一样
t.window：控制趋势项变化的速度,较小的值意味着更快的变化速度

data(nottem)
head(nottem,24)

## [1] 40.6 40.8 44.4 46.7 54.1 58.5 57.7 56.4 54.3 50.5 42.9 39.8 44.2 39.8 45.1
## [16] 47.0 54.1 58.7 66.3 59.9 57.0 54.2 39.7 42.8

stl(nottem, s.window = "periodic", t.window = 50)

## Call:
## stl(x = nottem, s.window = "periodic", t.window = 50)
##
## Components
## seasonal trend remainder
## Jan 1920 -9.3542441 49.94158 0.01266149
## Feb 1920 -9.8589768 49.89992 0.75905251
## Mar 1920 -6.8533560 49.85827 1.39509007
## Apr 1920 -2.7528761 49.81661 -0.36373152
## May 1920 3.5211867 49.77495 0.80386403
## Jun 1920 8.9998830 49.73329 -0.23317393
## Jul 1920 12.8601452 49.68456 -4.84470197
## Aug 1920 11.4833902 49.63582 -4.71921287
## Sep 1920 7.4481630 49.58709 -2.73525153
## Oct 1920 0.4621026 49.53835 0.49954308