抱歉,您的浏览器无法访问本站

本页面需要浏览器支持(启用)JavaScript


了解详情 >

HermitLSR

[曾记否,到中流击水,浪遏飞舟?]

  学习《时间序列分析及应用·Jonathan D.Cryer & Kung-Sik Chan 著》一书时的相关笔记——对于确定性趋势的参数估计,大致有五种:常均值模型、线性模型、二次式模型、季节均值模型、余弦模型。虽然每种模型各有特点和要求,但对于参数求解,一般都是使用OLSE。其区别仅在与设计矩阵的构造。本文将介绍这五种模型参数估计的R语言自编程序的实现。

常均值模型

  1. 参数估计的R语言函数

  常均值模型为:$$Y_t=\mu+X_t$$ 其中对所有的 ttE(X)=0E(X)=0。使用观测到的时间序列 Y1,Y2,,YnY_1,Y_2,…,Y_n 来估计μ\mu 。则可由OLSEOLSE得出μ\mu的估计为样本均值或平均数,定义为$$\bar{Y}=\frac{1}{n}\sum_{t=1}^n Y_t$$在模型假设成立的前提下,可以看出E(Yˉ)=μE(\bar{Y})=\mu,基于此常均值模型的参数估计算法函数ConReg()如下:

1
2
3
4
5
6
ConReg <- function(Y){

mu = sum(Y) / length(Y)
out = list(Y=Y, mu=mu)

}
  1. 参数估计的R语言函数检验

  由于课本并未给出相关数据,下面人为生成一组数据进行验证,做法如下:令模型中的μ=2\mu=2XtN(0,1)X_t\sim N(0,1),通过模型生成一组Yt,t=1,,100Y_t,t=1,…,100,以此使用上述函数ConReg()估计μ\mu的值,比较真实的μ=2\mu=2,判断算法的正确性。代码如下:

1
2
3
4
5
6
set.seed(1)
mu = 2
X <- rnorm(100,0,1)
Y <- mu + X + rnorm(100, 0, .2)
muest <- ConReg(Y)
muest$mu

  可以计算得估计μ=2.101\mu=2.101,与真实值相差不大,可以认为函数ConReg()是合理的。

线性模型

  1. 参数估计的R语言函数

  对于线性模型$$\mu_t=\beta_0+\beta_1t+\beta_2t+…+\beta_nt$$ 类似的可通过OLSEOLSE,计算得起参数估计β\beta应该为:$$\beta=(X’X)^{-1}XY$$基于此给出线性模型参数估计的模型LinReg()代码如下:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
LinReg <- function(X, y, lambda = 0, intercept = T){
n = dim(X)[1]
p = dim(X)[2]

#Coefficients
if(intercept == T){
X. = cbind(rep(1,n), X)
}else{
X. = X
}

beta = solve(t(X.)%*%X. + lambda*diag(dim(X.)[2]), t(X.)%*%y)

if (intercept == T){
beta1 = beta[2:(p+1)]
beta0 = beta[1]
}else{
beta1 = beta
beta0 =NULL
}

out = list(X=X, y=y, beta=beta1, beta0 = beta0,
intercept = intercept )
}

  1. 参数估计的R语言函数检验

  类似于常均值模型,人为产生数据进行验证,代码如下:

1
2
3
4
5
6
7
8
9
10
11
12
n <- 100
p <- 6
set.seed(1)
X <- matrix(rnorm(n*p), n, p)
beta = c(3, 1, rep(0,p-2))
beta0 = 0
set.seed(1)
y = X%*%beta + beta0 + rnorm(100, 0, .2)

result <- LinReg(X, y)
result$beta
result$beta0

  可以看出估计的参数β^\hat\beta与真实值非常接近,认为该函数没有问题。

二次模型

  1. 参数估计的R语言函数

  对于二次模型我们可以将其转化为线性模型进行参数估计,一般的二次模型如下:$$\mu_t=\beta_0+\beta_1t^2$$我们可以令X0=1,X1=t2X_0=1,X_1=t^2,这样将其转为,设计矩阵为X=(X0X1)X=(\begin{matrix}X_0&X_1 \end{matrix})的线性模型,为方便以后的函数调用,在设计R语言函数时推广为多项式模型,代码如下:

1
2
3
4
5
6
7
8
9
10
11
12
PolyDes <- function(t){
n <- length(t[, 1])
p <- length(t[1, ])
X <- matrix(nrow = n, ncol = p)
for (i in 1:p) {
X[, i] <- t[, i]^i
}

out = X

}

  1. 参数估计的R语言函数检验

  函数PolyDes()仅给出多项式模型的设计矩阵,要求其参数估计可再次调用线性模型函数LinRes()。以下给出模型的验证:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
beta <- c(5, 4, 3, 2, 1)
beta0 <- 0

t <- matrix(nrow = 100, ncol = 5)
set.seed(1)
r <- matrix(rnorm(100*5), 100, 5)
for (i in 1:5) {
t[, i] = 6-i + r[, i]
}

y = beta0 + t[,1]*beta[1] + t[,2]^2*beta[2] + t[,3]^3*beta[3] + t[,4]^4*beta[4] + t[,5]^5*beta[5]+rnorm(100, 0, .2)

X <- PolyDes(t)

result <- LinReg(X,y)
result$beta
result$beta0

  对比估计的参数和真实的参数β\beta是很接近的,可以认为函数是有效的。

季节趋势模型

  1. 参数估计的R语言函数

  对季节性趋势模型,其观测到的数据可表示为$$Y_t=\mu_t+X_t$$,其中,对所有t,E(Xt)=0E(X_t)=0。对于月度季节新数据,对μt\mu_t最常用的假设是存在12个常数(参数)β1,β2,,β12\beta_1,\beta_2,\cdots,\beta_{12},它们给出了12个月的期望平均气温。可以记为:

μ={β1t=1,13,15,β2t=2,14,16,β12t=12,24,36,\mu=\begin{cases}\beta_1&t=1,13,15,…\\ \beta_2&t=2,14,16,…\\ \vdots \\ \beta_{12}&t=12,24,36,…\end{cases}

该模型有时被称为季节均值模型。
  对于季节均值模型的估计,关键点在于将模型按月进行参数估计,这需要对设计矩阵XX进行特殊处理。这里的处理是,设计矩阵XX是大小为n×12n×12的示性矩阵,n为季节性时间序列的长度。XX的每一行就为一个示性函数,仅有唯一一个位置为1,即对于确定的iiXijX_{ij}表示第ii个数据是第jj月的。值得注意的是对于有截距项的模型如下:

μt=β0+β2I(t%12=2)++β12I(t%12=0)\mu_t=\beta_0+\beta_2I(t\%12=2)+\cdots+\beta_{12}I(t\%12=0)

参考前面的线性参数进估计,季节均值模型参数估计代码如下:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
SeaReg <- function(y, intercept = F){
if (!is.ts(y))
stop("y need to be a time series")
n <- length(y)
X <- matrix(ncol = 12, nrow = n)

month <-c("January", "February", "March", "April", "May", "June",
"July", "August", "September", "October", "November", "December")
month.y <- as.character(season(y))
for (i in 1:n) {
for (j in 1:12) {
if(month.y[i]==month[j]){
X[i,j] = 1
}else{
X[i,j] = 0
}
}
}

if(intercept == F){
X. = X
}else{
X[, 1]=rep(1, n)
X. = X

}

beta = solve(t(X.)%*%X., t(X.)%*%y)
out = list(y=y, X=X., beta=beta,intercept = intercept )
}

  1. 参数估计的R语言函数检验

  为对季节性趋势模型函数SeaReg()进行检验,采用TSA包中的tempdub数据进行检验。代码如下:

1
2
3
4
5
6
7
8
library(TSA)
data("tempdub")
#无截距项
result1 <- SeaReg(tempdub)
result1$beta
#有截距项
result2 <- SeaReg(tempdub, intercept = T)
result2$beta

  对比《时间序列》P24页的数据可知,有截距项和无截距项的参数估计的结果是基本一致的,故而该函数是可以接受的。

余弦趋势模型

  1. 参数估计的R语言函数

  余弦趋势模型的最简单形式可表示为:$$\mu_t=\beta_0+\beta_1cos(2\pi ft)+\beta_2sin(2\pi ft)$$类比二项式模型可知,其可也转换为线性模型求解,只不过需要对其设计矩阵做一定的设置,代码如下:

1
2
3
4
5
6
7
8
9
10
11
12
13
CosDes <- function(y, f=1){
if (!is.ts(y))
stop("y need to be a time series")
n = length(y)
t = as.numeric(time(y))
X1 = cos(2*pi*f*t)
X2 = sin(2*pi*f*t)

X = cbind(X1, X2)
out = X

}

  1. 参数估计的R语言函数检验

  采用TSA包中的tempdub数据进行检验。使用默认频率f=1,代码如下:

1
2
3
4
5
6
library(TSA)
data("tempdub")
X <- CosDes(tempdub, 1)
result<-LinReg(X,tempdub)
result$beta
result$beta0

  可以看到其参数估计的结果与课本P25页表3-5的数据一致,说明该函数的确可以完成余弦趋势模型的设计矩阵变换。




萌ICP备20201324号 | 渝ICP备2021006569号 | 公安备50011702500564号

又拍云logo     腾讯云logo