3

I'm attempting to run a bayesian time series analysis using the bsts package, but I'm getting bizarre forecasts when I set family to "poisson." My time series are counts and follow a Poisson distribution, so this is necessary. I suspect I need to establish the state specifications differently, but information is sparse online.

Example data:

as.integer(data)

[1]  NA  NA  NA  NA  NA  NA  NA  NA  NA   0  NA   0   1   0   0   0   6   3   7   6   2   0   1  20   3  14 326 702 500 47 157

Code used:

### Creating state specifications

ss=AddLocalLinearTrend(list(),y=data)
ss=AddSeasonal(ss,data,nseasons = 31)

### Model
model=bsts(data,ss,niter=500,family="poisson",na.action=na.omit)

### Predicting next value
pred=predict(model, horizon = 1)
pred

Output:

$`mean`
[1] 121474737

$median
[1] 0

$interval
        [,1]
2.5%           0
97.5% 2113254279

$distribution
         [,1]
[1,]          0
[2,]          0
[3,]          0
[4,]          0
[5,]  386555543
[6,]          0
[7,] 1315215054
[8,]          0
[9,]          0
[10,]          0
[11,]          0
[12,]          0
[13,]          0
[14,]          0
[15,]          0
[16,]          0
[17,]          0
[18,]          0
[19,]          0
[20,]          0
[21,]          0
[22,]          0
[23,]          0
[24,]          0
[25,]          0
[26,]          0
[27,]          0
[28,]          0
[29,]          0
[30,]        378
[31,]          0
[32,]          0
[33,]          0
[34,]          0
[35,]          0
[36,]          0
[37,]          0
[38,]          0
[39,] 1655342080
[40,]          0
[41,]          0
[42,]          0
[43,]          0
[44,]          0
[45,]          0
[46,]          0
[47,]          0
[48,]          0
[49,]          0
[50,]          0
[51,]          0
[52,] 1452117468
...

Any suggestions would be greatly appreciated! Thanks.

sebp
  • 1,787
  • 13
  • 24
user218440
  • 31
  • 4

1 Answers1

3

This is not well documented, but the trick is to take the log of your outcome data when adding things like trends and seasonality. You shouldn't take the log when passing in the data to the bsts function though because it will give you an error about not having integer inputs. For your example this would look like this:


### Creating state specifications

ss=AddLocalLinearTrend(list(),y=log1p(data))
ss=AddSeasonal(ss,log1p(data),nseasons = 31)

### Model
model=bsts(data,ss,niter=500,family="poisson",na.action=na.omit)

### Predicting next value
pred=predict(model, horizon = 1)
pred

Note: The log1p is important here to handle situations with 0 events in a given time period, it's just taking the log of data + 1 rather than data.

With your input data the output should look something like this (won't be exactly the same as the seed isn't set):

$mean
[1] 16651.53

$median
[1] 378.5

$interval
          [,1]
2.5%       2.0
97.5% 154583.7

$distribution
          [,1]
  [1,]    7090
  [2,]     282
  [3,]     125
  [4,]      96
  [5,]    1296
  [6,]      84
  [7,]       0
  [8,]     101
  [9,]     523
 [10,]      23
 [11,]    2357
 [12,]       8
 [13,]     917
 [14,]     309
 [15,]      50
 [16,]     181
 [17,]     134
 [18,]    9949
 [19,]     158
 [20,]      77
 [21,]     480
 [22,]    1936
 [23,]     475
 [24,]      59
 [25,]      11
 [26,]      90
 [27,]     162
 [28,]       4
 [29,]      44
 [30,]     150
 [31,]    1838
 [32,]    1720
 [33,]    1313
 [34,]      26
 [35,]      70
 [36,]     210
 [37,]    3667
 [38,]      66
 [39,]     259
 [40,]     111
 [41,]     727
 [42,]  178638
 [43,]    1329
 [44,]   86233
 [45,]     509
 [46,]    1696
 [47,]       2
 [48,]    9122
 [49,]      57
 [50,]     543
 [51,]     289
 [52,]    1355
 [53,]      12
 [54,]      86
 [55,]    1750
 [56,]     503
 [57,]   48816
 [58,]    1575
 [59,]    2942
 [60,]   16534
 [61,]   65270
 [62,]     742
 [63,]     105
 [64,]     194
 [65,]  797984
 [66,]     134
 [67,]     669
 [68,]     144
 [69,]    1013
 [70,]     199
 [71,]      46
 [72,]    3228
 [73,]   16962
 [74,]     221
 [75,]     475
 [76,]    2919
 [77,]     159
 [78,]      52
 [79,]      75
 [80,]      56
 [81,]       9
 [82,]    1680
 [83,]     282
 [84,]    1772
 [85,]    5569
 [86,]      43
 [87,]      23
 [88,]     292
 [89,]      13
 [90,]    1786
 [91,]     170
 [92,]     953
 [93,]      44
 [94,]     439
 [95,]      12
 [96,]     111
 [97,]      95
 [98,]     382
 [99,]     408
[100,]       4
[101,]    1426
[102,]      88
[103,]     189
[104,]    4735
[105,]      18
[106,]      32
[107,]       4
[108,]      77
[109,]    3760
[110,]      21
[111,]      36
[112,]    2628
[113,]     195
[114,]     385
[115,]     534
[116,]   34162
[117,]      36
[118,]     281
[119,]    7725
[120,]    6126
[121,]       9
[122,]     869
[123,]    4326
[124,]    5664
[125,]       6
[126,]    3354
[127,]      42
[128,]     483
[129,]    2709
[130,]     270
[131,]      25
[132,]    1422
[133,]    7915
[134,]     348
[135,]    2672
[136,]     385
[137,]     106
[138,]   55652
[139,]       3
[140,]     194
[141,]     106
[142,]       8
[143,]    3345
[144,]    5259
[145,]     533
[146,]       5
[147,]   24033
[148,]     255
[149,]    8830
[150,]     384
[151,]       8
[152,]    1466
[153,]  407763
[154,]     217
[155,]    1648
[156,]    2331
[157,]     157
[158,]    9887
[159,]   42085
[160,]     468
[161,]    1873
[162,]   10824
[163,]     329
[164,]    3091
[165,]      19
[166,]     749
[167,]      85
[168,]    1029
[169,]    1889
[170,]   20991
[171,]       1
[172,]      45
[173,]    1851
[174,]     103
[175,]    1331
[176,]    5674
[177,]     265
[178,]     102
[179,]       5
[180,]   14724
[181,]    1434
[182,]     131
[183,]    4588
[184,]       3
[185,]      31
[186,]      11
[187,]    1561
[188,]      38
[189,]   56054
[190,]   21187
[191,]       2
[192,]       2
[193,]   13736
[194,]      54
[195,]      41
[196,]    6497
[197,]    1738
[198,]     520
[199,]     575
[200,]     454
[201,]    2225
[202,]  118035
[203,]  122359
[204,]     959
[205,]   43948
[206,]     385
[207,]    3484
[208,]    3749
[209,]    1900
[210,]     224
[211,]     401
[212,]     543
[213,]     301
[214,]    1618
[215,]      30
[216,]     109
[217,]     165
[218,]     411
[219,]     211
[220,]    4775
[221,]      20
[222,]     359
[223,]      27
[224,]      57
[225,]      54
[226,]    1013
[227,]  135821
[228,]  184834
[229,]       3
[230,]     151
[231,]     130
[232,]    3008
[233,]    2281
[234,]    4501
[235,]      23
[236,]       0
[237,]     433
[238,]  106177
[239,]     191
[240,]    2019
[241,]       6
[242,]     502
[243,]      19
[244,]       8
[245,]    9798
[246,]      15
[247,]     879
[248,]     375
[249,]      76
[250,]     168
[251,]       2
[252,]  899456
[253,]      95
[254,]     311
[255,] 1158136
[256,]    1390
[257,]     281
[258,]      25
[259,]    1722
[260,]       1
[261,]   17394
[262,]    2789
[263,]       7
[264,]   18077
[265,]     117
[266,]   32022
[267,]     470
[268,]       0
[269,]     374
[270,]       0
[271,]       0
[272,]       2
[273,]     364
[274,]     346
[275,]     317
[276,]      54
[277,]       1
[278,]  160031
[279,]      84
[280,]  447847
[281,]     400
[282,]      23
[283,]      32
[284,]      14
[285,]     246
[286,]   95959
[287,]    1584
[288,]     107
[289,]      95
[290,]    4307
[291,]     555
[292,]    5822
[293,]       6
[294,]    2095
[295,]   25783
[296,]     843
[297,]   25186
[298,]     413
[299,]    6045
[300,]     190
[301,]    1854
[302,]    2408
[303,]     244
[304,]   11229
[305,]    1688
[306,]      16
[307,]   54187
[308,]    1276
[309,]     741
[310,]      41
[311,]   14750
[312,]       3
[313,]    1055
[314,]   19068
[315,]      89
[316,]    3524
[317,]     669
[318,]    2736
[319,]   35246
[320,]     166
[321,]     670
[322,]      20
[323,]      10
[324,]   10030
[325,]      81
[326,]    3416
[327,]     314
[328,]     199
[329,]     359
[330,]     197
[331,]   21780
[332,]     111
[333,]  258354
[334,]       6
[335,]    1773
[336,]    7888
[337,]  315910
[338,]      58
[339,]    9974
[340,]     313
[341,]     840
[342,]    3264
[343,]     280
[344,]     163
[345,]      14
[346,]    1926
[347,]     635
[348,]     177
[349,]     426
[350,]    8808
[351,]      17
[352,]     753
[353,]      20
[354,]      68
[355,]     600
[356,]     415
[357,]  162939
[358,]    3187
[359,]      24
[360,]     524
[361,]      44
[362,]   12529
[363,]    1102
[364,]   12555
[365,]     569
[366,]     240
[367,]     189
[368,]    1923
[369,]       6
[370,]      19
[371,]     191
[372,]    4633
[373,]      28
[374,]    1447
[375,]    1535
[376,]    4648
[377,]   40958
[378,]   62949
[379,]  212563
[380,]    2935
[381,]       1
[382,]   19464
[383,]   12656
[384,]      95
[385,]       7
[386,]    1247
[387,]   16725
[388,]   33275
[389,]      13
[390,]      69
[391,]     230
[392,]      38
[393,]      33
[394,]   54081
[395,]     218
[396,]     147
[397,]       7
[398,]      25
[399,]      27
[400,]    1271
[401,]       6
[402,]      47
[403,]   87372
[404,]     943
[405,]     851
[406,]      28
[407,]       8
[408,]   10685
[409,]     165
[410,]      21
[411,]    1383
[412,]    2847
[413,]     117
[414,]    2776
[415,]       0
[416,]     106
[417,]    6530
[418,]      48
[419,]     390
[420,]      74
[421,]     245
[422,]       1
[423,]      72
[424,]      67
[425,]    8450
[426,]     111
[427,]   83061
[428,]    2375
[429,]      39
[430,]    4657
[431,]     132
[432,]       6
[433,]     280
[434,]     159
[435,]       4
[436,]     540
[437,]      12
[438,]    3092
[439,]     447
[440,]     784
[441,]    2069
[442,]     241
[443,]      13
[444,]    7041
[445,]      51
[446,]       2
[447,]      13
[448,]     505
[449,]     618
[450,]      41

$original.series
  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 
 NA  NA  NA  NA  NA  NA  NA  NA  NA   0  NA   0   1   0   0   0   6   3   7   6   2   0   1  20   3  14 326 702 500 
 30  31 
 47 157 

attr(,"class")
[1] "bsts.prediction"
```
Pbrooks
  • 31
  • 2