1
$\begingroup$

I am trying to forecast a time series using the bsts package. The data is as follows:

enter image description here

The data is collected monthly and the goal is to forecast y by v1, v2, v3, v4, v6 and v7. y is a count variable, so I use poisson family.

library(bsts)
library(dplyr)
data_f<-structure(list(v1 = c(17242647.0588235, 18498333.3333333, 21890000, 
                          27110416.6666667, 34084800, 39341904.7619048, 45212500, 35353684.2105263, 
                          35353684.2105263, 35811304.3478261, 40493913.0434783, 44562500, 
                          46071176.4705882, 48077600, 45101904.7619048, 42889230.7692308, 
                          40886400, 40281250, 39787916.6666667, 39927391.3043478, 43662800, 
                          47327307.6923077, 50170833.3333333, 58340476.1904762, 62048421.0526316, 
                          66089583.3333333, 72500000, 92497307.6923077, 104050400, 112614000, 
                          137931739.130435, 125291304.347826, 113013076.923077, 111041880, 
                          108333800, 105900409.090909, 106789000, 99218400, 105324800, 
                          106098000, 113881000, 119995000), v2 = c(5.5e+07, 5.9e+07, 6.5e+07, 
                                                                   6.9e+07, 7.4e+07, 8.1e+07, 8.6e+07, 9.2e+07, 9.5e+07, 9.8e+07, 
                                                                   9.9e+07, 1.09e+08, 1.13e+08, 1.27e+08, 1.33e+08, 1.33e+08, 1.3e+08, 
                                                                   1.26e+08, 1.27e+08, 1.25e+08, 1.35e+08, 1.39e+08, 1.43e+08, 1.57e+08, 
                                                                   1.53e+08, 1.67e+08, 1.91e+08, 2.01e+08, 2.31e+08, 2.43e+08, 2.67e+08, 
                                                                   2.72e+08, 2.69e+08, 2.74e+08, 2.84e+08, 3.02e+08, 2.93e+08, 2.87e+08, 
                                                                   2.96e+08, 3.05e+08, 3.09e+08, 3.17e+08), v3 = c(1582, 1582, 1582, 
                                                                                                                   1672, 1672, 1672, 1764, 1764, 1764, 1882, 1882, 1882, 1979, 1979, 
                                                                                                                   1979, 2126, 2126, 2126, 2156, 2211, 2262, 2472, 2472, 2472, 3467, 
                                                                                                                   3467, 3467, 3467, 3467, 3467, 3467, 3467, 3467, 3467, 3467, 3467, 
                                                                                                                   3489, 3722, 3835, 3956, 3985, 4001), v4 = c(51663, 60976, 66582, 
                                                                                                                                                               81250, 102708, 120493, 121996, 123500, 111061, 108454, 118300, 
                                                                                                                                                               130323, 134135, 144001, 133362, 125505, 118750, 113765, 113928, 
                                                                                                                                                               114579, 126722, 130753, 135349, 150653, 155444, 160430, 174200, 
                                                                                                                                                               209365, 215418, 237377, 292610, 268402, 255041, 250225, 238348, 
                                                                                                                                                               243862, 243910, 221980, 234970, 247850, 259710, 274700), v6 = c(0.081, 
                                                                                                                                                                                                                               0.08, 0.082, 0.087, 0.097, 0.113, 0.134, 0.156, 0.18, 0.206, 
                                                                                                                                                                                                                               0.235, 0.269, 0.306, 0.342, 0.376, 0.404, 0.422, 0.427, 0.42, 
                                                                                                                                                                                                                               0.411, 0.4, 0.386, 0.37, 0.348, 0.322, 0.298, 0.278, 0.264, 0.258, 
                                                                                                                                                                                                                               0.26, 0.272, 0.29, 0.305, 0.322, 0.342, 0.364, 0.495, 0.469, 
                                                                                                                                                                                                                               0.476, 0.428, 0.432, 0.437), v7 = c(10000, 10000, 10000, 10000, 
                                                                                                                                                                                                                                                                   10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 
                                                                                                                                                                                                                                                                   10000, 10000, 10000, 10000, 10000, 10000, 10000, 30000, 30000, 
                                                                                                                                                                                                                                                                   30000, 30000, 30000, 30000, 30000, 30000, 30000, 30000, 30000, 
                                                                                                                                                                                                                                                                   30000, 30000, 30000, 30000, 30000, 30000, 30000, 30000, 30000, 
                                                                                                                                                                                                                                                                   30000, 30000), y = c(10, 20, 39, 21, 18, 10, 30, 28, 19, 15, 
                                                                                                                                                                                                                                                                                        16, 19, 7, 15, 10, 14, 11, 11, 14, 52, 8, 10, 8, 11, 7, 14, 11, 
                                                                                                                                                                                                                                                                                        15, 14, 15, 13, 12, 13, 17, 22, 18, 0, 0, 0, 0, 0, 0)), row.names = c(NA, 
                                                                                                                                                                                                                                                                                                                                                              42L), class = "data.frame")
data_f_train<-data_f[1:36,]
data_f_test<-data_f[37:42,]
data_f_test<-data_f_test %>% select(-y)
ss=AddLocalLinearTrend(list(),y=log1p(data_f_train$y))
ss=AddSeasonal(ss,log1p(data_f_train$y),nseasons = 12, season.duration = 1)
model2=bsts(y~.
        ,state.specification =ss,niter=100,family="poisson",data=data_f_train)

pred <- predict(model2,  newdata=data_f_test)
par(mfrow=c(2,1))
plot(pred,ylim=c(0,100))
data_f_test2<-data_f_test+10000
pred2 <- predict(model2,  newdata=data_f_test2)
plot(pred2,ylim=c(0,100))

It seems my steps is wrong completely since the forecast does not change as I change the covariate variables.

enter image description here

> head(model2$coefficients)
     (Intercept) v1 v2 v3 v4 v6 v7
[1,]           0  0  0  0  0  0  0
[2,]           0  0  0  0  0  0  0
[3,]           0  0  0  0  0  0  0
[4,]           0  0  0  0  0  0  0
[5,]           0  0  0  0  0  0  0
[6,]           0  0  0  0  0  0  0
Improve this question
$\endgroup$
1
  • $\begingroup$ Is this because the y variable is small but the covariate variables are too big?! $\endgroup$
    – Masoud
    Nov 3, 2021 at 16:45

1 Answer 1

Reset to default
2
+50
$\begingroup$

You are fitting a poisson model for count data, which means the predicted value can't be $< 0$. If there is a negative relationship between a covariate and the outcome, and the predicted count is already $0$, increasing the value of the covariate won't change the predicted value, since it can't go any lower.

To test this, you could try fitting the same data with a Gaussian linear model, probably by setting family = 'gaussian', or similar. This allows negative predicted values, and so the predictions should no longer be constant for these cases.

Improve this answer
$\endgroup$
4
  • $\begingroup$ I checked it, but the result is same( the result not change when covariates changes) $\endgroup$
    – Masoud
    Nov 15, 2021 at 19:16
  • $\begingroup$ Hi Masoud. I've just looked at the actual code you've posted. The y column, which you've annotated as "forecast" is actually the data, not the model predictions, so I'm afraid I don't understand your question, sorry. $\endgroup$
    – Eoin
    Nov 15, 2021 at 23:01
  • $\begingroup$ We want to predict 'y'(row 37-42 in dataframe data_f) , based on 'V_i'. $\endgroup$
    – Masoud
    Nov 19, 2021 at 16:51
  • $\begingroup$ Could you explain? 'The y column, which you've annotated as "forecast" is actually the data, not the model predictions'. $\endgroup$
    – Masoud
    Nov 22, 2021 at 5:27

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.