更新于 

确定性趋势参数估计

常均值模型

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

  常均值模型为：

Yt=μ+XtY_t=\mu+X_t

其中对所有的tt有E(X)=0E(X)=0。使用观测到的时间序列Y1,Y2,…,YnY_1,Y_2,…,Y_n 来估计μ\mu 。则可由OLSEOLSE得出μ\mu的估计为样本均值或平均数，定义为

Yˉ=1n∑t=1nYt\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=2，Xt∼N(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语言函数

  对于线性模型

μt=β0+β1t+β2t+…+βnt\mu_t=\beta_0+\beta_1t+\beta_2t+…+\beta_nt

类似的可通过OLSEOLSE，计算得起参数估计β\beta应该为：

β=(X′X)−1XY\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语言函数

  对于二次模型我们可以将其转化为线性模型进行参数估计，一般的二次模型如下：

μt=β0+β1t2\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，即对于确定的ii，XijX_{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语言函数

  余弦趋势模型的最简单形式可表示为:

μt=β0+β1cos(2πft)+β2sin(2πft)\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的数据一致，说明该函数的确可以完成余弦趋势模型的设计矩阵变换。