确定性趋势参数估计
常均值模型
参数估计的R语言函数
常均值模型为:
Y t = μ + X t Y_t=\mu+X_t
Y t = μ + X t
其中对所有的t t t 有E ( X ) = 0 E(X)=0 E ( X ) = 0 。使用观测到的时间序列Y 1 , Y 2 , … , Y n Y_1,Y_2,…,Y_n Y 1 , Y 2 , … , Y n 来估计μ \mu μ 。则可由O L S E OLSE O L S E 得出μ \mu μ 的估计为样本均值或平均数,定义为
Y ˉ = 1 n ∑ t = 1 n Y t \bar{Y}=\frac{1}{n}\sum_{t=1}^n Y_t
Y ˉ = n 1 t = 1 ∑ n Y t
在模型假设成立的前提下,可以看出E ( Y ˉ ) = μ E(\bar{Y})=\mu E ( Y ˉ ) = μ ,基于此常均值模型的参数估计算法函数ConReg()
如下:
1 2 3 4 5 6 ConReg <- function ( Y) { mu = sum ( Y) / length ( Y) out = list ( Y= Y, mu= mu) }
参数估计的R语言函数检验
由于课本并未给出相关数据,下面人为生成一组数据进行验证,做法如下:令模型中的μ = 2 \mu=2 μ = 2 ,X t ∼ N ( 0 , 1 ) X_t\sim N(0,1) X t ∼ N ( 0 , 1 ) ,通过模型生成一组Y t , t = 1 , … , 100 Y_t,t=1,…,100 Y t , t = 1 , … , 1 0 0 ,以此使用上述函数ConReg()
估计μ \mu μ 的值,比较真实的μ = 2 \mu=2 μ = 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 μ = 2 . 1 0 1 ,与真实值相差不大,可以认为函数ConReg()
是合理的。
线性模型
参数估计的R语言函数
对于线性模型
μ t = β 0 + β 1 t + β 2 t + … + β n t \mu_t=\beta_0+\beta_1t+\beta_2t+…+\beta_nt
μ t = β 0 + β 1 t + β 2 t + … + β n t
类似的可通过O L S E OLSE O L S E ,计算得起参数估计β \beta β 应该为:
β = ( X ′ X ) − 1 X Y \beta=(X'X)^{-1}XY
β = ( X ′ X ) − 1 X Y
基于此给出线性模型参数估计的模型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 ] 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 ) }
参数估计的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 β ^ 与真实值非常接近,认为该函数没有问题。
二次模型
参数估计的R语言函数
对于二次模型我们可以将其转化为线性模型进行参数估计,一般的二次模型如下:
μ t = β 0 + β 1 t 2 \mu_t=\beta_0+\beta_1t^2
μ t = β 0 + β 1 t 2
我们可以令X 0 = 1 , X 1 = t 2 X_0=1,X_1=t^2 X 0 = 1 , X 1 = t 2 ,这样将其转为,设计矩阵为X = ( X 0 X 1 ) X=(\begin{matrix}X_0&X_1 \end{matrix}) X = ( X 0 X 1 ) 的线性模型,为方便以后的函数调用,在设计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 }
参数估计的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 β 是很接近的,可以认为函数是有效的。
季节趋势模型
参数估计的R语言函数
对季节性趋势模型,其观测到的数据可表示为$$Y_t=\mu_t+X_t$$,其中,对所有t,E ( X t ) = 0 E(X_t)=0 E ( X t ) = 0 。对于月度季节新数据,对μ t \mu_t μ t 最常用的假设是存在12个常数(参数)β 1 , β 2 , ⋯ , β 12 \beta_1,\beta_2,\cdots,\beta_{12} β 1 , β 2 , ⋯ , β 1 2 ,它们给出了12个月的期望平均气温。可以记为:
μ = { β 1 t = 1 , 13 , 15 , … β 2 t = 2 , 14 , 16 , … ⋮ β 12 t = 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}
μ = ⎩ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎧ β 1 β 2 ⋮ β 1 2 t = 1 , 1 3 , 1 5 , … t = 2 , 1 4 , 1 6 , … t = 1 2 , 2 4 , 3 6 , …
该模型有时被称为季节均值 模型。
对于季节均值模型的估计,关键点在于将模型按月进行参数估计,这需要对设计矩阵X X X 进行特殊处理。这里的处理是,设计矩阵X X X 是大小为n × 12 n×12 n × 1 2 的示性矩阵,n为季节性时间序列的长度。X X X 的每一行就为一个示性函数,仅有唯一一个位置为1,即对于确定的i i i ,X i j X_{ij} X i j 表示第i i i 个数据是第j j j 月的。值得注意的是对于有截距项的模型如下:
μ t = β 0 + β 2 I ( t % 12 = 2 ) + ⋯ + β 12 I ( t % 12 = 0 ) \mu_t=\beta_0+\beta_2I(t\%12=2)+\cdots+\beta_{12}I(t\%12=0)
μ t = β 0 + β 2 I ( t % 1 2 = 2 ) + ⋯ + β 1 2 I ( t % 1 2 = 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 ) }
参数估计的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页的数据可知,有截距项和无截距项的参数估计的结果是基本一致的,故而该函数是可以接受的。
余弦趋势模型
参数估计的R语言函数
余弦趋势模型的最简单形式可表示为:
μ t = β 0 + β 1 c o s ( 2 π f t ) + β 2 s i n ( 2 π f t ) \mu_t=\beta_0+\beta_1cos(2\pi ft)+\beta_2sin(2\pi ft)
μ t = β 0 + β 1 c o s ( 2 π f t ) + β 2 s i n ( 2 π f t )
类比二项式模型可知,其可也转换为线性模型求解,只不过需要对其设计矩阵做一定的设置,代码如下:
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 }
参数估计的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的数据一致,说明该函数的确可以完成余弦趋势模型的设计矩阵变换。