Прогноз зависимых значений с mvrnorm и включают временную автокорреляцию - PullRequest
Новая
       28

Прогноз зависимых значений с mvrnorm и включают временную автокорреляцию

0 голосов
/ 04 июля 2018

У меня есть матрицы со значениями (вес, зрелость и т. Д.) По временному шагу и возрастному классу, и я хотел бы прогнозировать будущие значения неопределенно. Возрастные классы не являются независимыми, поэтому я использовал 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

Хотя последняя часть кажется не совсем правильной.

1 Ответ

0 голосов
/ 04 июля 2018

Первая часть этого ответа - как объяснить временную автокорреляцию в исходном вопросе. Вторая часть добавляет ответ о многовариантном случае для пересмотренного вопроса.

Часть 1: 

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 / веб / пакеты / CARBayesST / виньетки / CARBayesST.pdf

...