我有按时间步长和年龄类划分的值(权重、成熟度等(的矩阵，我想不确定地预测未来的值。年龄等级不是独立的，所以我一直在使用 mvrnorm 来处理这个问题。如何在预测中也获得(滞后 1(时间自相关？
这是我想在R中做的事情：
library(MASS)
# dummy matrix: rows are time steps columns are dependent classes (ages)
x <- matrix(rnorm(20),4,5,dimnames = list(years=c('year1','year2','year3','year4'),ages=c('age1','age2','age3','age4','age5')))
# what I have so far to get next year's values (the goal would be to predict several years)
sigma <- cov(x) #covariance matrix
delta <- mvrnorm(1,rep(0,ncol(x)),cov(x)) # deviations
xl <- tail(x,1) #last year values
xp <- xl+delta #new values
# There is no temporal autocorrelation in here of course
xnew <- rbind(x,xp)
matplot(xnew,type='l')
# So I would need new values based on something like this:
rho <- apply(x,2,function(x) acf(x)$acf[2,1,1])
delta <- mvrnorm(1,xl,cov(x))
xp <- rho*xl+(1-rho)*delta
不过最后一部分感觉不对。
这个答案的第一部分是如何解释原始问题中的时间自相关。第 2 部分根据修订后的问题添加了有关多变量案例的答案。
篇一：
library(MASS)
# dummy matrix: rows are time steps columns are dependent classes (ages)
x <- matrix(rnorm(20),4,5)
# what I have so far to get next year's values (the goal would be to predict several years)
sigma <- cov(x)
delta <- mvrnorm(1,rep(0,ncol(x)),cov(x))
xl <- tail(x,1)
xp <- xl+delta #new values
# There is no temporal autocorrelation in here of course
xnew <- rbind(x,xp)
matplot(xnew,type='l')
# Clean up / construct your data set
dat <- as.data.frame(x)
dat$year <- c(2014,2015,2016,2017)
dat <- rbind(dat, c(xp, 2018))
colnames(dat) <- c("maturity", "age", "height", "sales", "year")
# Account for Temporal Autocorrelation
library(nlme)
mdl.ac <- gls(sales ~ year, data=dat,
correlation = corAR1(form=~year),
na.action=na.omit)
summary(mdl.ac)
Generalized least squares fit by REML Model: sales ~ year Data: dat AIC BIC logLik 14.01155 10.406 -3.005773 Correlation Structure: ARMA(1,0) Formula: ~year Parameter estimate(s): Phi1 0.1186508 Coefficients: Value Std.Error t-value p-value (Intercept) 1.178018 0.5130766 2.2959883 0.1054 year 0.012666 0.3537748 0.0358023 0.9737 Correlation: (Intr) year 0.646 Standardized residuals: 1 2 3 4 5 0.3932124 -0.4053291 -1.8081473 0.0699103 0.8821300 attr(,"std") [1] 0.5251018 0.5251018 0.5251018 0.5251018 0.5251018 attr(,"label") [1] "Standardized residuals" Residual standard error: 0.5251018 Degrees of freedom: 5 total; 3 residual
第 2 部分：
# Account for Temporal Autocorrelation
library(nlme)
mdl.ac <- gls(year ~ height + sales + I(maturity*age), data=dat,
correlation = corAR1(form=~year),
na.action=na.omit)
summary(mdl.ac)
Generalized least squares fit by REML Model: year ~ height + sales + I(maturity * age) Data: dat AIC BIC logLik 15.42011 3.420114 -1.710057 Correlation Structure: ARMA(1,0) Formula: ~year Parameter estimate(s): Phi1 0 Coefficients: Value Std.Error t-value p-value (Intercept) 0.2100381 0.4532345 0.4634203 0.7237 height -0.7602539 0.7758925 -0.9798444 0.5065 sales -0.1840694 0.8327382 -0.2210411 0.8615 I(maturity * age) 0.0449278 0.1839260 0.2442712 0.8475 Correlation: (Intr) height sales height -0.423 sales 0.214 -0.825 I(maturity * age) 0.349 -0.941 0.889 Standardized residuals: 1 2 3 4 5 -7.004956e-17 -4.985525e-01 -1.319137e+00 -1.568271e+00 -1.441708e+00 attr(,"std") [1] 0.3962277 0.3962277 0.3962277 0.3962277 0.3962277 attr(,"label") [1] "Standardized residuals" Residual standard error: 0.3962277 Degrees of freedom: 5 total; 1 residual
另请参阅CARBayesST及其小插图以获取替代方法：
https://cran.r-project.org/web/packages/CARBayesST/vignettes/CARBayesST.pdf