Ensemble with bootstrap forecasts for a transformed variable

[ac_1]

Hi Tom,

Questions regarding BACK_TRANSFORMED, mean, median and SE in an ensemble with bootstrap forecasts for the LOG of a variable ly.

I can generate BOOTSTRAP forecasts for ly:

Code: Select all

declare vector[series] fcast_D(ndraws*2)
comp Nfcast_D = %rows(fcast_D)
declare vector[series] efcast_D(ndraws*2)
comp Nefcast_D = %rows(efcast_D)

do draw = 1, ndraws
   boxjenk(noprint,constant,ar=||1||,diffs=1,ma=0,define=EnsembleEq,MAXL) ly istart iend
   uforecast(equation=EnsembleEq,BOOTSTRAP) fore_D iend+1 iend+h
   set fcast_D(draw) iend+1 iend+h = fore_D
   set efcast_D(draw) iend+1 iend+h = exp(fore_D); * back-transformed median
end do draw
do draw = ndraws+1, (ndraws*2)
   boxjenk(noprint,constant,ar=||1,2,3,4,5||,diffs=1,ma=0,define=EnsembleEq,MAXL) ly istart iend
   uforecast(equation=EnsembleEq,BOOTSTRAP) fore_D iend+1 iend+h
   set fcast_D(draw) iend+1 iend+h = fore_D
   set efcast_D(draw) iend+1 iend+h = exp(fore_D); * back-transformed median
end do draw

prin / fcast_D efcast_D

And PI's for the BACK-TRANSFORMED BOOTSTRAP's:

Code: Select all

set efcastMIN_D iend+1 iend+h = work=%xt(efcast_D,t),%fractiles(work,||.00||)(1); * 100%

....

set efcastmean_D iend+1 iend+h = work=%xt(efcast_D,t),%avg(work); * mean
set efcastmedian_D iend+1 iend+h = work=%xt(efcast_D,t),%fractiles(work,||.50||)(1); * median

....

set efcastMAX_D iend+1 iend+h = work=%xt(efcast_D,t),%fractiles(work,||1.00||)(1); * 100%

efcastmean_D appears to be the correct BACK-TRANSFORMED BIAS-ADJ MEAN.


The following calculates the MEAN's and SE's for the TRANSFORMED series ly:

Code: Select all

set acsum iend+1 iend+h = 0.0
set acsum_SS iend+1 iend+h = 0.0
do draw = 1, (ndraws*2)
   set acsum iend+1 iend+h = acsum+fcast_D(draw)
   set acsum_SS iend+1 iend+h = acsum_SS+(fcast_D(draw))^2

* save fcast_D(draw) in a cloud of simulations
   dec vect[series] columns_D(ndraws*2)
   set columns_D(draw) iend+1 iend+h = fcast_D(draw)
end do draw

prin / fcast_D columns_D; * same

set res_D_M iend+1 iend+h = acsum/(ndraws*2)
set res_D_MSS iend+1 iend+h = acsum_SS/(ndraws*2)
set res_D_SE iend+1 iend+h = sqrt(res_D_MSS-(res_D_M)^2)

prin / res_D_M res_D_SE

Using: res_D_M and res_D_SE,
how do I calculate
(a) the same BACK-TRANSFORMED BIAS-ADJ MEAN: efcastmean_D?
(b) and the correct (aggregated) SE's: res_D_SE, for the TRANSFORMED series ly?

as

Code: Select all

set fmean_D iend+1 iend+h = exp(res_D_M) + 0.5*(res_D_SE**2.0)*exp(res_D_M)

is not the same as efcastmean_D?


thanks,
Amarjit

[TomDoan]

1. The mean of a log normal is exp(mu+.5sigma^2) where mu and sigma are for the underlying Normal. That's not what you're computing.
2. That relationship between the moments of X and the moments of exp(X) is specific to X being Normal, and is only true in population, not in samples, and not with arbitrary distributions.
3. Why are you aggregating statistics from two different models? 1 through ndraw is one model and ndraw+1 to 2*ndraw is another.

[ac_1]

Are the forecasts being generated correctly as a cloud of bootstraps?

1. The mean of a log normal is exp(mu+.5sigma^2) where mu and sigma are for the underlying Normal. That's not what you're computing.

Code: Select all

set fmean_D iend+1 iend+h = exp(res_D_M) + 0.5*(res_D_SE**2.0)*exp(res_D_M)

is the calculation for the BACK-TRANSFORMED BIAS-ADJ MEAN,

yhat = exp(forecast logged series) + 0.5*(stderrs squared logged series)*exp(forecast logged series) i.e. Taylor series expansion of a back-transformed log transformation.


Also, is efcastmean_D the right BACK-TRANSFORMED BIAS-ADJ MEAN in the PI's calculation?

2. That relationship between the moments of X and the moments of exp(X) is specific to X being Normal, and is only true in population, not in samples, and not with arbitrary distributions.

Yes.

3. Why are you aggregating statistics from two different models? 1 through ndraw is one model and ndraw+1 to 2*ndraw is another.

I am producing a bootstrap ensemble model. 50% of the 1st and 50% of the 2nd. And then attempting to aggregate the statistics as above.

Where in the acsum calculations (which I think are correct for a single series, but may not be for an ensemble as I would be aggregating SE's from two different models) are the errors?

[TomDoan]

1. The mean of a log normal is exp(mu+.5sigma^2) where mu and sigma are for the underlying Normal. That's not what you're computing.

Code: Select all

set fmean_D iend+1 iend+h = exp(res_D_M) + 0.5*(res_D_SE**2.0)*exp(res_D_M)

is the calculation for the BACK-TRANSFORMED BIAS-ADJ MEAN,

yhat = exp(forecast logged series) + 0.5*(stderrs squared logged series)*exp(forecast logged series) i.e. Taylor series expansion of a back-transformed log transformation.

I'm not sure what the point of that is. If you're getting your forecasts by taking the mean of a simulated/bootstrapped set of samples, why would the (at best) approximately biased-transformed estimate of the forecasts of the series in the "wrong" form be of any interest?

Also, is efcastmean_D the right BACK-TRANSFORMED BIAS-ADJ MEAN in the PI's calculation?

I'm not sure I would describe it that way (there is no "bias-adjustment"). But it is a bootstrapped estimate of the forecast of the levels.

3. Why are you aggregating statistics from two different models? 1 through ndraw is one model and ndraw+1 to 2*ndraw is another.

I am producing a bootstrap ensemble model. 50% of the 1st and 50% of the 2nd. And then attempting to aggregate the statistics as above.

Where in the acsum calculations (which I think are correct for a single series, but may not be for an ensemble as I would be aggregating SE's from two different models) are the errors?

Throwing them together in a calculation presupposes that equal weights are correct. But there is no difference between that and (equally) weighting the separately estimated forecasts, while doing separate calculations allows you to easily change the relative weights. And aggregating them for computing 2nd moments is just wrong.

[ac_1]

3. Why are you aggregating statistics from two different models? 1 through ndraw is one model and ndraw+1 to 2*ndraw is another.

I am producing a bootstrap ensemble model. 50% of the 1st and 50% of the 2nd. And then attempting to aggregate the statistics as above.

Where in the acsum calculations (which I think are correct for a single series, but may not be for an ensemble as I would be aggregating SE's from two different models) are the errors?

Throwing them together in a calculation presupposes that equal weights are correct. But there is no difference between that and (equally) weighting the separately estimated forecasts, while doing separate calculations allows you to easily change the relative weights. And aggregating them for computing 2nd moments is just wrong.

Combination forecasts is in Enders(2014) AETS 4thEdn p.109. IMHO a simple average is reasonable as optimal weights may be an "overfit".

How do I calculate the forecast mean's and forecast SE's (aggregating them) from the acsum calculations?

[TomDoan]

I'm not sure that figuring out weights from just a pair of models has any chance of being an "overfit", and if one of the two models really isn't very good, you are going to get trash forecasts with no possibility of fixing it. Now there is no practical way to compute truly "optimal" weights for a sizable collection of forecast procedures because you can't reasonably calculate a full covariance matrix that would be needed to minimize the variance of a linear combination, but a 2 x 2 matrix wouldn't be a stretch. However, there are methods for weighting forecasts based upon the behavior of each forecast procedure in isolation so you only need K statistics rather than K(K+1)/2.

Assuming, however, that you want to do equal weights, your forecasts are 1/2(f1-hat+f2-hat). If you are getting f1-hat and f2-hat by averaging across bootstrapped draws, you can pair the draws and compute the sample average of squared errors from the average of an f1 simulation and an f2 simulation. You *cannot* compute any useful 2nd moment statistic about of a big string of 2 x draws bootstrap values---even though the overall forecast will be the average across the 2 x draws bootstrap values, the second moment calculations on separate f1 and f2 errors can't be aggregated.

[ac_1]

Thanks!

I'm not sure that figuring out weights from just a pair of models has any chance of being an "overfit", and if one of the two models really isn't very good, you are going to get trash forecasts with no possibility of fixing it. Now there is no practical way to compute truly "optimal" weights for a sizable collection of forecast procedures because you can't reasonably calculate a full covariance matrix that would be needed to minimize the variance of a linear combination, but a 2 x 2 matrix wouldn't be a stretch. However, there are methods for weighting forecasts based upon the behavior of each forecast procedure in isolation so you only need K statistics rather than K(K+1)/2.

What are those methods?

Assuming, however, that you want to do equal weights, your forecasts are 1/2(f1-hat+f2-hat). If you are getting f1-hat and f2-hat by averaging across bootstrapped draws, you can pair the draws and compute the sample average of squared errors from the average of an f1 simulation and an f2 simulation. You *cannot* compute any useful 2nd moment statistic about of a big string of 2 x draws bootstrap values---even though the overall forecast will be the average across the 2 x draws bootstrap values, the second moment calculations on separate f1 and f2 errors can't be aggregated.

For model 1 and model 2, to calculate
log forecasts: f_1_hat & f_2_hat
back-transformed median forecasts: ef_1_hat & ef_2_hat

I would:

Code: Select all

do draw = 1, ndraws
  boxjenk(noprint,constant,ar=||1||,diffs=1,ma=0,define=EEq_1,MAXL) ly istart iend
  uforecast(equation=EEq_1,BOOTSTRAP) fore_1 iend+1 iend+h
  set fcast_1(draw) iend+1 iend+h = fore_1
  set efcast_1(draw) iend+1 iend+h = exp(fore_1)

  boxjenk(noprint,constant,ar=||1,2,3,4,5||,diffs=1,ma=0,define=EEq_2,MAXL) ly istart iend
  uforecast(equation=EEq_2,BOOTSTRAP) fore_2 iend+1 iend+h
  set fcast_2(draw) iend+1 iend+h = fore_2
  set efcast_2(draw) iend+1 iend+h = exp(fore_2)
end do draw


do draw = 1, ndraws
   set sum_1 iend+1 iend+h = sum_1+fcast_1(draw)
   set esum_1 iend+1 iend+h = esum_1+fcast_1(draw)

   set sum_2 iend+1 iend+h = sum_2+fcast_2(draw)
   set esum_2 iend+1 iend+h = esum_2+fcast_2(draw)
end do draw

set f_1_hat iend+1 iend+h = sum_1/ndraws
set ef_1_hat iend+1 iend+h = esum_1/ndraws
set f_2_hat iend+1 iend+h = sum_2/ndraws
set ef_2_hat iend+1 iend+h = esum_2/ndraws

Correct?


Then how to "pair the draws and compute the sample average of squared errors from the average of an f1 simulation and an f2 simulation" in RATS?

Typically, SE's are required for PI's and to plot fan charts; and can be calculated for a single series (no (ndraws*2)) using the acsum calculations. Right?

There must be a way to generate fan charts from ensembling/combining model 1 and model 2 after averaging across bootstrapped draws - maybe using?

Code: Select all

set fcast80L_D iend+1 iend+h = work=%xt((f_1_hat+f_2_hat)/2.0,t),%fractiles(work,||.1||)(1)
set fcast80U_D iend+1 iend+h = work=%xt((f_1_hat+f_2_hat)/2.0,t),%fractiles(work,||.9||)(1)

[TomDoan]

Thanks!

I'm not sure that figuring out weights from just a pair of models has any chance of being an "overfit", and if one of the two models really isn't very good, you are going to get trash forecasts with no possibility of fixing it. Now there is no practical way to compute truly "optimal" weights for a sizable collection of forecast procedures because you can't reasonably calculate a full covariance matrix that would be needed to minimize the variance of a linear combination, but a 2 x 2 matrix wouldn't be a stretch. However, there are methods for weighting forecasts based upon the behavior of each forecast procedure in isolation so you only need K statistics rather than K(K+1)/2.

What are those methods?

See, e.g. Zou, Hui & Yang, Yuhong, 2004. "Combining time series models for forecasting," International Journal of Forecasting, Elsevier, vol. 20(1), pages 69-84. I'm sure you can find others.

Assuming, however, that you want to do equal weights, your forecasts are 1/2(f1-hat+f2-hat). If you are getting f1-hat and f2-hat by averaging across bootstrapped draws, you can pair the draws and compute the sample average of squared errors from the average of an f1 simulation and an f2 simulation. You *cannot* compute any useful 2nd moment statistic about of a big string of 2 x draws bootstrap values---even though the overall forecast will be the average across the 2 x draws bootstrap values, the second moment calculations on separate f1 and f2 errors can't be aggregated.

For model 1 and model 2, to calculate
log forecasts: f_1_hat & f_2_hat

back-transformed median forecasts: ef_1_hat & ef_2_hat

I would:

Code: Select all

do draw = 1, ndraws
  boxjenk(noprint,constant,ar=||1||,diffs=1,ma=0,define=EEq_1,MAXL) ly istart iend
  uforecast(equation=EEq_1,BOOTSTRAP) fore_1 iend+1 iend+h
  set fcast_1(draw) iend+1 iend+h = fore_1
  set efcast_1(draw) iend+1 iend+h = exp(fore_1)

  boxjenk(noprint,constant,ar=||1,2,3,4,5||,diffs=1,ma=0,define=EEq_2,MAXL) ly istart iend
  uforecast(equation=EEq_2,BOOTSTRAP) fore_2 iend+1 iend+h
  set fcast_2(draw) iend+1 iend+h = fore_2
  set efcast_2(draw) iend+1 iend+h = exp(fore_2)
end do draw


do draw = 1, ndraws
   set sum_1 iend+1 iend+h = sum_1+fcast_1(draw)
   set esum_1 iend+1 iend+h = esum_1+fcast_1(draw)

   set sum_2 iend+1 iend+h = sum_2+fcast_2(draw)
   set esum_2 iend+1 iend+h = esum_2+fcast_2(draw)
end do draw

set f_1_hat iend+1 iend+h = sum_1/ndraws
set ef_1_hat iend+1 iend+h = esum_1/ndraws
set f_2_hat iend+1 iend+h = sum_2/ndraws
set ef_2_hat iend+1 iend+h = esum_2/ndraws

Correct?


Then how to "pair the draws and compute the sample average of squared errors from the average of an f1 simulation and an f2 simulation" in RATS?

Typically, SE's are required for PI's and to plot fan charts; and can be calculated for a single series (no (ndraws*2)) using the acsum calculations. Right?

There must be a way to generate fan charts from ensembling/combining model 1 and model 2 after averaging across bootstrapped draws - maybe using?

Code: Select all

set fcast80L_D iend+1 iend+h = work=%xt((f_1_hat+f_2_hat)/2.0,t),%fractiles(work,||.1||)(1)
set fcast80U_D iend+1 iend+h = work=%xt((f_1_hat+f_2_hat)/2.0,t),%fractiles(work,||.9||)(1)

A bootstrapped draw for the combined model is just .5*(efcast_1(draw)+efcast_2(draw)). Save that when you're saving the other two.

[ac_1]

Thanks. Yes - this works.

I would like to plot the distribution at EACH time step for a Mixture distribution from BOOTSTRAP forecasts. Would that be using SCATTER and DENSITY as?

density(smoothing=1.5) fcasts_D(1) 1 ndraws xf1 f1_D
scatter(key=upright,head="", footer="Bootstrap Estimation Density Forecasts/Conditional PDFs fcasts_D(1)") 1
# xf1 f1_D / 1
etc etc

And if the forecast means are not from BOOTSTRAP forecasts, then either the:
- analytical normal PDF's
- simulating 1000000 random normal variates PDF's
as per here viewtopic.php?f=33&t=3638, as I would expect "smoothed" forecast distributions.

OR is there a way to plot a Mixture Normal distribution?


Regarding 2nd moment statistic/SE, how about

Var(aX + bY) = a^2*Var(X) + b^2*Var(Y) + 2*a*b*cor(X,Y)*sqrt(Var(X))*sqrt(Var(Y))

for an ENSEMBLE?

[TomDoan]

An average is not a "mixture". A mixture of Normals takes a probabilistic selection among separate Normal distributions. That's very, very different from what you are doing.

Why would you simulate a Normal distribution? All you would be doing is injecting simulation error.

It is true that var(A+B)=var(A)+2cov(A,B)+var(B). But that has no value unless you can compute cov(A,B), which you can't with the information you have. Again, you have a bootstrapped draw for the average by taking the average of a pair of bootstraps for the separate models---you don't have to re-invent anything.

[ac_1]

Another idea for aggregating SE's: https://en.wikipedia.org/wiki/Inverse-v ... _weighting

[TomDoan]

Another idea for aggregating SE's: https://en.wikipedia.org/wiki/Inverse-v ... _weighting

Aggregating SE's?? I'm not sure what you mean by that. At any rate, that doesn't really apply to what you are doing since it assumes the different measurements are independent. (That doesn't mean that it is an unreasonable way to feasibly weight forecasts, but the underlying theory doesn't apply). Note that this result is basically a special case of the weight calculation in Enders, but with the assumption of independence.

[ac_1]

There are 3 things here:

1. Optimal weighted forecasts

Enders has 4 methods:
(a) simple averaging i.e. equal weights
(b) weighting by the forecast error variances, eqn 2.74
(c) the regression method, eqn 2.75
(d) constructing the weights using SBC, p.111

The RATS example file enders4p111.rpf; results

Code: Select all

(b) eqn 2.74
             ar7   ar6   ar2  ar127 ar1ma1 ar2ma1 ar2ma17
var(eit)    0.186 0.180 0.170 0.170  0.169  0.177   0.169
varInv(eit) 5.366 5.544 5.878 5.878  5.918  5.639   5.931
wi          0.134 0.138 0.146 0.146  0.147  0.140   0.148

(c) eqn 2.75
Constrained-regression based weights 0.00000 0.00000 0.00000 0.32597 0.17046 0.00000 0.50357

(a) is straightforward.
(d) I haven't constructed the SBC weights, yet.


2. Mixture distribution: Previously, I was trying to generate a mixture distribution (1 to ndraws) and (ndraws+1 to ndraws*2) i.e. 50% model 1and 50% model 2.

3. Aggregating SE's: There must be a way?

[TomDoan]

Re (2). No. It's not a mixture distribution. Mixture has a specific meaning which I've explained. Equal averaging is (a).

(3). Compute using the training sample. It's a single methodology which can be applied to any method of computing and combining forecasts. And no. In general there is no other way. Standard errors of a linear combination require the covariances which in general don't exist from the estimation of the models.

[ac_1]

Re (2). No. It's not a mixture distribution. Mixture has a specific meaning which I've explained. Equal averaging is (a).

(3). Compute using the training sample. It's a single methodology which can be applied to any method of computing and combining forecasts. And no. In general there is no other way. Standard errors of a linear combination require the covariances which in general don't exist from the estimation of the models.

Thanks. I can do (1), and (3) Var(aX+bY) for a Normal, haven't attempted a Log-Normal or a Non-Central Chi-Squared.

The question is how to generate a mixture x% from model A and (1-x%) from model B: and forecast (i) analytically and (ii) bootstrapped?