I am trying to implement boostrap prediction interval example of FPP3 book in python for learning purposes (https://otexts.com/fpp3/prediction-intervals.html).
Prediction interval is estimated by simulating 5000 sample paths along forecasting horizon (h = 10). For each horizon and each path, random error (== residual in this case) value is sampled from collection of 252 known values and added to the previous forecasted value, leading to a random path each time. After that, for each horizon, lower and upper percentiles of the values of 5000 paths are calculated to get prediction interval at certain probability (i.e. 95%).
Data set: 2015 Google stock closing price, 252 businness days. Forecast horizon: 10 days. Model: Naive. Residuals: 252 values -> (observation value - naive model prediction) for each day.
bootstrap =  for i in range(0, 5000): path =  for j in range(0, 10): e = random.choice(list(residuals.values)) if (j == 0): path[j] = path[j] + e else: path[j] = path[j-1] + e bootstrap.append(path.to_numpy()) df = pd.DataFrame(bootstrap)
Finally, I have dataframe 'df':
Prediction intervals are calculated based on percentiles for each forecasting horizon:
path_l =  path_h =  ser = df.quantile(0.025, axis = 0) ser = ser.set_axis(path_l.index) path_l.update(ser) ser = df.quantile(0.975, axis = 0) ser = ser.set_axis(path_h.index) path_h.update(ser)
Normality check of distribution of simulated values (paths) for each forecasting horizon:
for i in range(0, 10): a = df[i].tolist() figure = plt.figure() figure.suptitle('histogram_bootstrap_step_' + str(i+1)) plt.hist(a, bins = 30) print(str(i+1) + ". step") print(stats.jarque_bera(a)) print(anderson(a))
Regarding bootstrapped prediction intervals, according to the histogram of samples drawn of residuals at a certain forecasting horizon (lets say h = 10 , last forecast), it seems that values could be normally distributed. However, normality tests (Anderson-Darling, Jarques-Bera) definitely show that distribution of drawn residuals is not normal, statistic values are much higher than required to reject the null hypothesis:
10. step Jarque_beraResult(statistic=1390.292955526538, pvalue=0.0) AndersonResult(statistic=47.28770856219671, critical_values=array([0.576, 0.655, 0.786, 0.917, 1.091]), significance_level=array([15. , 10. , 5. , 2.5, 1. ]))
Idk if the non-normality causes problem when taking percentiles as lower/upper prediction intervals.
Shortly, question is that is it required the bootrapped prediction interval to be normally distributed? Last graph ("Bootstrap intervals from the naïve method") on page:
Thank you very much.