Добавление отстающих переменных в модель lm?
Я использую lm для временного ряда, который работает довольно хорошо на самом деле, и он супер супер быстро.
Скажем, моя модель:
> formula <- y ~ x
Я тренирую это на тренировочном наборе:
> train <- data.frame( x = seq(1,3), y = c(2,1,4) )
> model <- lm( formula, train )
... и я могу сделать прогнозы для новых данных:
> test <- data.frame( x = seq(4,6) )
> test$y <- predict( model, newdata = test )
> test
x y
1 4 4.333333
2 5 5.333333
3 6 6.333333
Это прекрасно работает, и это очень быстро.
Я хочу добавить в модель отстающие переменные. Теперь я мог бы сделать это, увеличив свой первоначальный набор упражнений:
> train$y_1 <- c(0,train$y[1:nrow(train)-1])
> train
x y y_1
1 1 2 0
2 2 1 2
3 3 4 1
обновить формулу:
formula <- y ~ x * y_1
... и тренировка будет работать нормально:
> model <- lm( formula, train )
> # no errors here
Однако проблема заключается в том, что не существует способа использовать "предсказывать", потому что нет способа заполнения y_1 в тестовом наборе пакетным образом.
Теперь, для многих других регрессионных вещей, есть очень удобные способы выразить их в формуле, например poly(x,2) и т.д., и они работают напрямую, используя немодифицированные учебные и тестовые данные.
Итак, мне интересно, есть ли способ выражения отстающих переменных в формуле, чтобы можно было использовать predict? В идеале:
formula <- y ~ x * lag(y,-1)
model <- lm( formula, train )
test$y <- predict( model, newdata = test )
... без необходимости увеличивать (не уверены, правильно ли это слово) учебные и тестовые наборы данных и просто иметь возможность напрямую использовать predict?
Ответы
Ответ 1
Посмотрите, например. dynlm, который дает вам операторы с задержкой. В более общем плане в представлении задач по эконометрике и временным рядам будет намного больше возможностей для просмотра.
Вот начало его примеров - отставание в один и двенадцать месяцев:
R> data("UKDriverDeaths", package = "datasets")
R> uk <- log10(UKDriverDeaths)
R> dfm <- dynlm(uk ~ L(uk, 1) + L(uk, 12))
R> dfm
Time series regression with "ts" data:
Start = 1970(1), End = 1984(12)
Call:
dynlm(formula = uk ~ L(uk, 1) + L(uk, 12))
Coefficients:
(Intercept) L(uk, 1) L(uk, 12)
0.183 0.431 0.511
R>
Ответ 2
Следующее предложение Dirk на dynlm, я не мог понять, как предсказать, но поиск этого привел меня к пакету dyn через https://stats.stackexchange.com/questions/6758/1-step-ahead-predictions-with-dynlm-r-package
Затем, после нескольких часов экспериментов, я применил следующую функцию для обработки прогноза. На пути было довольно много "находок", например, вы не можете казаться временными рядами rbind, и результат предсказания компенсируется start и целым рядом таких вещей, поэтому я чувствую, что этот ответ добавляет значительно по сравнению с просто именованием пакета, хотя я получил ответ Дирка.
Итак, решение, которое работает, это:
- используйте пакет
dyn
- используйте следующий метод для прогнозирования
метод predDyn:
# pass in training data, test data,
# it will step through one by one
# need to give dependent var name, so that it can make this into a timeseries
predictDyn <- function( model, train, test, dependentvarname ) {
Ntrain <- nrow(train)
Ntest <- nrow(test)
# can't rbind ts apparently, so convert to numeric first
train[,dependentvarname] <- as.numeric(train[,dependentvarname])
test[,dependentvarname] <- as.numeric(test[,dependentvarname])
testtraindata <- rbind( train, test )
testtraindata[,dependentvarname] <- ts( as.numeric( testtraindata[,dependentvarname] ) )
for( i in 1:Ntest ) {
result <- predict(model,newdata=testtraindata,subset=1:(Ntrain+i-1))
testtraindata[Ntrain+i,dependentvarname] <- result[Ntrain + i + 1 - start(result)][1]
}
return( testtraindata[(Ntrain+1):(Ntrain + Ntest),] )
}
Пример использования:
library("dyn")
# size of training and test data
N <- 6
predictN <- 10
# create training data, which we can get exact fit on, so we can check the results easily
traindata <- c(1,2)
for( i in 3:N ) { traindata[i] <- 0.5 + 1.3 * traindata[i-2] + 1.7 * traindata[i-1] }
train <- data.frame( y = ts( traindata ), foo = 1)
# create testing data, bunch of NAs
test <- data.frame( y = ts( rep(NA,predictN) ), foo = 1)
# fit a model
model <- dyn$lm( y ~ lag(y,-1) + lag(y,-2), train )
# look at the model, it a perfect fit. Nice!
print(model)
test <- predictDyn( model, train, test, "y" )
print(test)
# nice plot
plot(test$y, type='l')
Вывод:
> model
Call:
lm(formula = dyn(y ~ lag(y, -1) + lag(y, -2)), data = train)
Coefficients:
(Intercept) lag(y, -1) lag(y, -2)
0.5 1.7 1.3
> test
y foo
7 143.2054 1
8 325.6810 1
9 740.3247 1
10 1682.4373 1
11 3823.0656 1
12 8686.8801 1
13 19738.1816 1
14 44848.3528 1
15 101902.3358 1
16 231537.3296 1
Изменить: хм, это очень медленно. Даже если я ограничу данные в subset постоянными несколькими строками набора данных, для прогнозирования требуется около 24 миллисекунд или для моей задачи 0.024*7*24*8*20*10/60/60= 1.792 hours: -O
Ответ 3
Вот мысль:
Почему бы вам не создать новый фрейм данных? Заполните кадр данных необходимыми регрессорами. У вас могут быть такие столбцы, как L1, L2,..., Lp для всех лагов любой переменной, которую вы хотите, и затем вы можете использовать свои функции точно так же, как и для регрессии поперечного сечения.
Потому что вам не придется работать с вашими данными каждый раз, когда вы вызываете функции фитинга и предсказания, но однажды преобразуете данные, это будет значительно быстрее. Я знаю, что Eviews и Stata обеспечивают отстающих операторов. Это правда, что есть какое-то удобство. Но это также неэффективно, если вам не нужны все функции, такие как "lm". Если у вас есть несколько сотен тысяч итераций для выполнения, и вам просто нужен прогноз или прогноз и значение информационных критериев, таких как BIC или AIC, вы можете победить "lm" в скорости, избегая вычислений, которые вы не будете use - просто напишите оценку OLS в функции, и вы хорошо пойдете.
