时间序列中的确定性趋势参数估计（R语言）

  对于确定性趋势的参数估计，大致有五种：常均值模型、线性模型、二次式模型、季节均值模型、余弦模型。虽然每种模型各有特点和要求，但对于参数求解，一般都是使用OLSE。其区别仅在与设计矩阵的构造。本文将介绍这五种模型参数估计的R语言自编程序的实现。
注：
1、本人最近正在学习《时间序列分析及应用》一书，本文的相关理论也来自于此书，水平有限，如有错误，还请多批评指正；
2、本文未贴出运行结果，有需要的朋友可以自行运行。（数据来源R包：TSA）

常均值模型

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

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

ConReg <- function(Y){

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

}

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

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

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

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

线性模型

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

  对于线性模型$${\mu }_t={\beta }_0+{\beta }_{1}t+{\beta }_{2}t+\dots +,{\beta }_{n}t$$ ,类似的可通过$OLSE$，计算得起参数估计$\beta$应该为：$$\beta =(X^{\prime }X{)}^{-1}XY$$基于此给出线性模型参数估计的模型LinReg()代码如下：

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 )
} 

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

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

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 }_1{t}^2$$我们可以令$X_0=1,X_1=t^2$，这样将其转为，设计矩阵为$X=(\begin{array}{cc}{X}_0& X_1\end{array})$的线性模型，为方便以后的函数调用，在设计R语言函数时推广为多项式模型，代码如下：

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
  
}

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

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

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(X_t)=0$。对于月度季节新数据，对${\mu }_t$最常用的假设是存在12个常数（参数）${\beta }_1,{\beta }_2,\cdots ,{\beta }_{12}$，它们给出了12个月的期望平均气温。可以记为：$$\mu =\{\begin{array}{ll}{\beta }_1& t=1,13,15,\dots \\ {\beta }_2& t=2,14,16,\dots \\ ⋮& \\ {\beta }_{12}& t=12,24,36,\dots \end{array}$$该模型有时被称为季节均值模型。
  对于季节均值模型的估计，关键点在于将模型按月进行参数估计，这需要对设计矩阵$X$进行特殊处理。这里的处理是，设计矩阵$X$是大小为$n\times 12$的示性矩阵，n为季节性时间序列的长度。$X$的每一行就为一个示性函数，仅有唯一一个位置为1，即对于确定的$i$，$X_{ij}$表示第$i$个数据是第$j$月的。值得注意的是对于有截距项的模型如下：$${\mu }_t={\beta }_0+{\beta }_{2}I(t\%12=2)+\cdots +{\beta }_{12}I(t\%12=0)$$参考前面的线性参数进估计，季节均值模型参数估计代码如下:

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 )
}

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

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

library(TSA)
data("tempdub")
#无截距项
result1 <- SeaReg(tempdub)
result1$beta
#有截距项
result2 <- SeaReg(tempdub, intercept = T)
result2$beta

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

余弦趋势模型

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

  余弦趋势模型的最简单形式可表示为:$${\mu }_t={\beta }_0+{\beta }_{1}cos(2\pi ft)+{\beta }_{2}sin(2\pi ft)$$类比二项式模型可知，其可也转换为线性模型求解，只不过需要对其设计矩阵做一定的设置，代码如下：

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

}

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

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

library(TSA)
data("tempdub")
X <- CosDes(tempdub, 1)
result<-LinReg(X,tempdub)
result$beta
result$beta0

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