Allin Cottrell and Ignacio Diaz-Emparanza 1.0 2026-11-07 Exponential smoothing: simple, Holt and Holt-Winters C53 ETS MAINWIN/Variable/Filter pdfdoc:ETS.pdf *ETS graphical interface* - The `model` selector chooses among simple exponential smoothing (level only), Holt's method (level and trend), and Holt-Winters with an additive or multiplicative seasonal component. The seasonal models require data with periodicity greater than 1 (quarterly, monthly, etc.). - The parameters α, β and γ are the smoothing coefficients for the level, trend and seasonal components of the time series, respectively. Only those relevant to the chosen model are used (α for simple smoothing, α and β for Holt, all three for Holt-Winters). Setting these values manually may be considered an "expert" (or pedagogical) option. If you leave the `parms` selection at `automatic` gretl will estimate them for you. - In the multiplicative Holt-Winters case the seasonal effect is proportional to the level-plus-trend component; in the additive case it's simply added to level-plus-trend. - The `ntest` selector allows you to reserve a certain number of observations (at the end of the available data) for testing forecasts. If the sample contains trailing observations beyond the range of the y data, setting this to zero means that out-of-sample forecasts will be made instead. - To see more details of the procedure, increase the `verbose` level. examples scalar fixed_theta = 0 scalar fixed_init = 0 string yname = argname(y, "y") matrix init = {} bundle args = ets_default_args() if exists(specs) fixed_theta = inbundle(specs, "theta") fixed_init = inbundle(specs, "init") # a seasonal model implies a trend unless told otherwise if inbundle(specs, "seasonal") && !inbundle(specs, "trend") args.trend = specs.seasonal > 0 endif errorif(bcheck(&args, specs), "aborting") endif ets_checks(args, fixed_theta, fixed_init) scalar td = args.trend scalar sf = args.seasonal scalar m = sf > 0 ? $pd : 0 scalar ntest = args.ntest matrix limits = ets_ranges(y, ntest, m > 0 ? 4 * m : 4) scalar t1 = limits[1] scalar t2 = limits[2] scalar t3 = limits[3] smpl t1 t3 scalar verbose = args.verbose series level = NA list L = y level if td series trend = NA L += trend endif if sf > 0 series seasonal = NA L += seasonal endif series forecast = NA L += forecast if verbose > 0 printf "%s\n", ets_model_name(td, sf) printf "Training range: %s to %s\n", obslabel(t1), obslabel(t2) if fixed_theta printf "Fixed parameters:\n" ets_print_theta(args.theta) endif if fixed_init printf "Fixed initializer:\n" ets_print_init(args.init, td, sf) endif if ntest > 0 printf "Testing forecast: %s to %s\n\n", obslabel(t2+1), obslabel(t3) elif t3 > t2 printf "Out-of-sample forecast: %s to %s\n\n", obslabel(t2+1), obslabel(t3) else printf "Number of forecasts: 0\n\n" endif endif if fixed_init init = vec(args.init) else init = ets_auto_init(y, td, sf) if verbose > 1 print init endif endif # training data as a column vector smpl t1 t2 matrix my = {y} smpl t1 t3 if fixed_theta matrix theta = vec(args.theta) if sf > 0 && !fixed_init err = ets_given_theta(theta, my, &init, td, sf, m, verbose) if err funcerr "Failed to optimize initializer" endif endif else matrix theta = ets_auto_parms(my, &init, td, sf, m, verbose, fixed_init) if !ok(theta[1]) funcerr "Failed to estimate parameters" endif endif theta = vec(theta) # training: states and one-step-ahead forecasts matrix Wtr = ets_filter(theta, init, my, td, sf, m) scalar fcol = 3 + td + (m > 0) matrix u = Wtr[,1] - Wtr[,fcol] scalar n = rows(u) scalar SSE = sumc(u.^2) scalar RMSE_train = sqrt(SSE / n) scalar MAE_train = sumc(abs(u)) / n # transcribe to series over t1..t2 smpl t1 t2 level = Wtr[,2] if td trend = Wtr[,3] endif if sf > 0 seasonal = Wtr[,3+td] endif forecast = Wtr[,fcol] # h-step-ahead forecasts over t2+1..t3 (test or out-of-sample) if t3 > t2 scalar H = t3 - t2 matrix hv = seq(1, H)' matrix fc = Wtr[n,2] * ones(H, 1) if td fc += hv * Wtr[n,3] endif if sf > 0 # last m seasonals, cycled out to H periods matrix sh = Wtr[n-m+1:n,3+td] * ones(1, ceil(H/m)) sh = vec(sh)[1:H] if sf == 1 fc += sh else fc = fc .* sh endif endif smpl t2+1 t3 forecast = fc endif smpl t1 t3 if verbose > 0 if verbose > 1 printf "\n" endif printf "training SSE = %g, RMSE = %g, MAE = %g\n", SSE, RMSE_train, MAE_train print L -o endif set skip_missing off matrix S = {L} nparms = ets_nparms(td, m, fixed_theta, fixed_init) bundle ret = _(theta, init, trend=td, seasonal=sf, ntest, S, RMSE_train, MAE_train, yname, limits, nparms) if ntest > 0 ets_add_eval_stats(&ret) if verbose > 0 printf "testing RMSE = %g, MAE = %g\n", ret.RMSE_test, ret.MAE_test endif endif if verbose > 0 printf "\n" endif if sf != 2 && t3 > t2 ets_add_fcse(&ret) endif if $version >= 20262 ets_add_virtual_series(&ret) endif return ret "simple exponential smoothing" "Holt (trend)" "Holt-Winters, additive" "Holt-Winters, multiplicative" "automatic" "pre-specified" level, α trend, β season, γ ntest, see Help string yname = argname(y, "y") if yname != "y" rename y @yname endif # decode the model selector trend = model > 1 seasonal = model > 2 ? model - 2 : 0 bundle ret if parms == 1 ret = ETS(@yname, _(trend, seasonal, ntest, verbose)) else matrix th3 = {alpha, beta, gamma} matrix theta = th3[1:1+trend+(seasonal>0)]' ret = ETS(@yname, _(theta, trend, seasonal, ntest, verbose)) endif return ret "actual and forecast" "components" list L = mat2list(b.S) if ptype == 1 string yname = b.yname if yname != "y" rename y @yname endif if b.limits[3] > b.limits[2] scalar c1 c2 hmax ets_get_ci_vals(b, ci_levels, &c1, &c2, &hmax) df = rows(b.S) - hmax - b.nparms if b.seasonal != 2 hbands = ets_ci_bands(df, c1, c2) else # simulation-based bands matrix Y = ets_mult_sim(&b, 10000, df) scalar T = rows(b.S) matrix B1 = ets_sim_band(Y, T, c1) matrix B2 = ets_sim_band(Y, T, c2) series ctr1 = B1[,1] series hw1 = B1[,2] series ctr2 = B2[,1] series hw2 = B2[,2] bundles hbands = array(2) hbands[1] = _(center="ctr2", width="hw2", style="fill", color="#BBD9F5") hbands[2] = _(center="ctr1", width="hw1", style="fill", color="#6AADEC") endif list P = @yname forecast plot P options --time-series --with-lines --bands=hbands printf "set title 'forecast with %d and %d percent intervals'", c1, c2 end plot --output=display else gnuplot @yname forecast --time --with-lines --output=display endif else series irreg = b.S[,1] - b.S[,ets_fcol(b)] string cs = "level" if b.trend cs ~= " trend" endif if b.seasonal > 0 cs ~= " seasonal" endif tsplots @cs irreg --output=display endif scalar hmax = b.limits[3] - b.limits[2] errorif(hmax <= 0, "No post-training forecasts are available") if exists(ci_levels) ets_check_ci_levels(ci_levels) else ci_levels = {80, 95} endif scalar df = rows(b.S) - hmax - b.nparms scalar ci1 = min(ci_levels) scalar ci2 = max(ci_levels) scalar T = rows(b.S) scalar fcol = ets_fcol(b) matrix fc = b.S[T-hmax+1:T,fcol] matrix CI = zeros(hmax, 4) if b.seasonal == 2 # simulation-based intervals matrix Y = ets_mult_sim(&b, 10000, df) matrix B1 = ets_sim_band(Y, T, ci1) matrix B2 = ets_sim_band(Y, T, ci2) CI[,1] = B1[T-hmax+1:T,1] - B1[T-hmax+1:T,2] CI[,2] = B1[T-hmax+1:T,1] + B1[T-hmax+1:T,2] CI[,3] = B2[T-hmax+1:T,1] - B2[T-hmax+1:T,2] CI[,4] = B2[T-hmax+1:T,1] + B2[T-hmax+1:T,2] else scalar f1 = critical(t, df, (100 - ci1) / 200) scalar f2 = critical(t, df, (100 - ci2) / 200) matrix me = b.S[T-hmax+1:T,fcol+1] CI[,1] = fc - f1 * me CI[,2] = fc + f1 * me CI[,3] = fc - f2 * me CI[,4] = fc + f2 * me endif strings cs = array(4) cs[1] = sprintf("Lo%d", ci1) cs[2] = sprintf("Hi%d", ci1) cs[3] = sprintf("Lo%d", ci2) cs[4] = sprintf("Hi%d", ci2) cnameset(CI, cs) return CI # number of estimated parameters: smoothing parameters plus free # initial-state elements (seasonals are constrained to sum to zero # or to average 1, hence m-1 of them). With fixed theta, a # non-seasonal initializer is not estimated either. scalar np = fixed_theta ? 0 : 1 + td + (m > 0) scalar ni = 1 + td + (m > 0 ? m - 1 : 0) if fixed_init || (fixed_theta && m == 0) ni = 0 endif return np + ni # column of the forecast in b.S return 3 + b.trend + (b.seasonal > 0) strings gk = defarray("α", "β", "γ") scalar n = nelem(theta) printf " " loop i=1..n printf "%s = %g", gk[i], theta[i] if i < n printf ", " endif endloop printf "\n" printf " level = %g", init[1] scalar i = 2 if td printf ", trend = %g", init[i++] endif printf "\n" if sf && $pd > 4 scalar n = nelem(init) printf " seasonals:\n" loop j=i..n printf " %+f\n", init[i++] endloop elif sf scalar n = nelem(init) printf " seasonals = {%g, ", init[i++] loop j=i..n printf "%g%s", init[i++], j < n ? ", " : "}" endloop printf "\n" endif printf "\n" matrix ret = zeros(1,3) catch smpl --contiguous y if $error funcerr "Missing values encountered" endif ret[1] = $t1 # start of contiguous y ret[2] = $t2 # end of contiguous y scalar tail = $tmax - $t2 if ntest > 0 ret[2] = $t2 - ntest ret[3] = $t2 elif tail > 0 ret[3] = $tmax else ret[3] = $t2 endif ntrain = ret[2] - ret[1] + 1 errorif(ntrain < minlen, "Insufficient observations for training") return ret scalar td = args.trend scalar sf = args.seasonal if sf > 0 errorif($pd == 1, "Periodicity = 1, no seasonality!") errorif(td == 0, "seasonal models require a trend") endif if fixed_theta scalar np = 1 + td + (sf > 0) errorif(nelem(args.theta) != np, sprintf("theta must be a %d-vector", np)) loop i=1..np errorif(args.theta[i] < 0 || args.theta[i] > 1, sprintf("theta[%d] is out of bounds", i)) endloop endif if fixed_init scalar ni = 1 + td + (sf > 0 ? $pd : 0) errorif(nelem(args.init) != ni, sprintf("init must be a %d-vector", ni)) errorif(sumc(ok(vec(args.init))) < ni, "init contains missing values") if sf == 2 matrix iv = vec(args.init) errorif(iv[1] <= 0 || minc(iv[3:ni]) <= 0, "multiplicative model: level and seasonals must be positive") endif endif if sf == 0 return td ? "Holt double exponential smoothing" : "Simple exponential smoothing" else return sprintf("Holt-Winters, using %s formula", sf == 1 ? "additive" : "multiplicative") endif Target series # Rule-of-thumb initial state, for period 0. Non-seasonal: initial # level is the first observation, initial trend the first # difference. Seasonal: first-period average, trend from the # first two major periods, and seasonal indices relative to the # first-period average. if sf == 0 matrix ret = zeros(1 + td, 1) ret[1] = y[$t1] if td ret[2] = y[$t1+1] - y[$t1] rnameset(ret, "level trend") else rnameset(ret, "level") endif return ret endif scalar m = $pd scalar t2 = $t1+m-1 scalar am1 = 0 scalar am2 = 0 matrix ret = zeros(m + 2, 1) # sums of y over the first two major periods loop t=$t1..t2 am1 += y[t] am2 += y[t + m] endloop ret[1] = am1 / m ret[2] = (am2 - am1) / (m * m) # seasonal indices i = 3 loop t=$t1..t2 if sf == 1 ret[i] = y[t] - ret[1] else ret[i] = y[t] / ret[1] endif i++ endloop strings S = array(2 + m) S[1] = "level" S[2] = "trend" loop i=1..m S[i+2] = sprintf("seas%d", i) endloop rnameset(ret, S) return ret scalar α = theta[1] scalar β = td ? theta[2] : 0 scalar γ = m > 0 ? theta[2+td] : 0 scalar n = 1 + td + m matrix A = zeros(n, n) A[1,1] = 1-α if td A[1,2] = 1-α A[2,1] = -α*β A[2,2] = 1-α*β endif if m > 0 A[1,n] = -α if td A[2,n] = -α*β endif scalar sr = 2 + td # row of the newest seasonal A[sr,1] = -γ*(1-α) if td A[sr,2] = -γ*(1-α) endif A[sr,n] = 1-γ*(1-α) loop j=sr+1..n # shift the older seasonals down one slot A[j,j-1] = 1 endloop endif return A # loading of y_t on the state (row vector) scalar α = theta[1] matrix d = {α} if td d ~= α * theta[2] endif if m > 0 d ~= theta[2+td] * (1-α) ~ zeros(1, m-1) endif return d # forecast selection vector matrix c = ones(1 + td, 1) if m > 0 c |= zeros(m-1, 1) | 1 endif return c # returns T x (3 + td + (m>0)) matrix # {y, level, [trend], [seasonal], forecast} scalar T = rows(my) scalar k = 1 + td # state order is newest seasonal first; @init is oldest first matrix z0 = init[1:k] if m > 0 z0 |= mreverse(init[k+1:k+m]) endif matrix Z = varsimul(ets_A(theta, td, m), my * ets_d(theta, td, m), z0') matrix ret = my ~ Z[2:T+1,1:k] if m > 0 ret ~= Z[2:T+1,k+1] endif ret ~= Z[1:T,] * ets_c(td, m) return ret # BFGS callback, additive case: SSE given @theta, with the initial # state concentrated out. Forecasts are linear in the initial # state, f = f0 + G z0, where f0 is the forecast path from a zero # state and row t of G is c'A^(t-1). Note that G can be had from # one varsimul() call with transition A' and initial value c. # # If seasonal, the seasonals are constrained to sum to zero: the # initial level and the seasonal mean are not separately identified # (adding k to the level and subtracting k from every seasonal # leaves all forecasts unchanged), so without this the regressors # would be collinear. scalar T = rows(my) scalar k = 1 + td scalar n = k + m matrix A = ets_A(theta, td, m) matrix c = ets_c(td, m) matrix Z = varsimul(A, my * ets_d(theta, td, m), zeros(1, n)) matrix f0 = Z[1:T,] * c matrix G = varsimul(A', zeros(T-1, n), c') matrix e, q if m > 0 # put seasonal columns in @init order (oldest first) matrix X = G[,1:k] ~ mreverse(G[,k+1:n]')' # impose seas_m = -(seas_1 + ... + seas_{m-1}) matrix Xs = X[,1:n-1] - X[,n] * (zeros(1, k) ~ ones(1, m-1)) q = mols(my - f0, Xs, &e) init[1:n-1] = q init[n] = -sumc(q[k+1:n-1]) else q = mols(my - f0, G, &e) init[1:k] = q endif return sumc(e.^2) # returns T x 5 matrix {y, level, trend, seasonal, forecast} scalar α = theta[1] scalar β = theta[2] scalar γ = theta[3] scalar m = nelem(init) - 2 scalar T = rows(my) # s[t] holds s_{t-m}; s[t+m] receives s_t matrix s = init[3:m+2] | zeros(T, 1) matrix ret = my ~ zeros(T, 4) scalar a_ = init[1] scalar b_ = init[2] scalar a = 0 scalar b = 0 scalar sm = 0 loop t=1..T sm = s[t] ret[t,5] = (a_ + b_) * sm a = α * my[t] / sm + (1-α) * (a_ + b_) b = β * (a - a_) + (1-β) * b_ s[t+m] = γ * my[t] / a + (1-γ) * sm ret[t,2] = a ret[t,3] = b a_ = a b_ = b endloop ret[,4] = s[m+1:m+T] return ret # map (α, β, γ, level, trend, seas_1..seas_{m-1}) to the full # initializer, with seasonals normalized to average 1: scaling # level and trend by k and the seasonals by 1/k leaves all # forecasts unchanged, so one seasonal is redundant matrix sv = par[6:m+4] return par[4:5] | sv | (m - sumc(sv)) # BFGS callback, multiplicative case matrix W = ets_filter_mult(par[1:3], ets_mult_init(par, m), my) return sumc((W[,1] - W[,5]).^2) if sf == 2 return ets_filter_mult(theta, init, my) else return ets_filter_add(theta, init, my, td, m) endif return ets_crit_add(theta, my, &init, td, m) / f0 # criterion (SSE) given both @theta and @init matrix W = ets_filter(theta, init, my, td, sf, m) return sumc((W[,1] - W[,cols(W)]).^2) return ets_crit_fixed(theta, init, my, td, sf, m) / f0 return ets_crit_mult(par, my, m) / f0 # Starting values for the smoothing parameters, by coarse grid # search over the criterion. The criterion can be quite irregular # (and the first L-BFGS-B step is a full gradient step, which for # a gradient of large magnitude can throw us onto a corner of the # box from which we cannot escape), so a fixed start such as # (0.5, 0.5) is not safe. On return @f0 holds the criterion value # at the chosen point, for use in scaling the criterion. scalar np = 1 + td + (m > 0) matrix P = {0.05; 0.1; 0.2; 0.4; 0.6; 0.8} if np > 1 matrix bg = {0.01; 0.05; 0.1; 0.3; 0.6} P = (P ** ones(rows(bg), 1)) ~ (ones(rows(P), 1) ** bg) endif if np > 2 matrix gg = {0.01; 0.1; 0.3} P = (P ** ones(rows(gg), 1)) ~ (ones(rows(P), 1) ** gg) endif scalar n = rows(P) matrix fv = zeros(n, 1) matrix tmp = init loop i=1..n if fixed_init fv[i] = ets_crit_fixed(P[i,]', init, my, td, sf, m) elif sf == 2 fv[i] = ets_crit_mult(P[i,]' | init[1:m+1], my, m) else fv[i] = ets_crit_add(P[i,]', my, &tmp, td, m) endif endloop scalar imin = iminr(fv') f0 = fv[imin] return P[imin,]' # driver for automatic estimation of parameters; if @fixed_init # is non-zero, @init is held constant and only theta is estimated scalar np = 1 + td + (m > 0) scalar f0 matrix b = ets_start_values(my, init, td, sf, m, fixed_init, &f0) if sf == 2 && !fixed_init # the additive cases solve for the initializer internally b |= init[1:m+1] endif scalar len = nelem(b) matrix bounds = zeros(len, 3) loop i=1..len if i <= np bounds[i,] = {i, 0, 1} else bounds[i,] = {i, -$huge, $huge} endif endloop if verbose == 3 set max_verbose full elif verbose == 2 set max_verbose on set bfgs_verbskip 5 endif if fixed_init catch scalar SSE = BFGScmin(&b, bounds, ets_crit_fixed_scaled(b, init, my, td, sf, m, f0)) elif sf == 2 catch scalar SSE = BFGScmin(&b, bounds, ets_crit_mult_scaled(b, my, m, f0)) else catch scalar SSE = BFGScmin(&b, bounds, ets_crit_add_scaled(b, my, &init, td, m, f0)) endif err = $error if err printf "Got error %d from BFGScmin\n", err b[1] = NA else # make sure @init corresponds to the final parameters if fixed_init # nothing to do elif sf == 2 init[1:m+2] = ets_mult_init(b, m) else SSE = ets_crit_add(b, my, &init, td, m) endif if verbose > 0 if verbose > 1 printf "\n" endif printf "estimated parameter values:\n" ets_print_theta(b[1:np]) if !fixed_init printf "\nestimated initializer:\n" ets_print_init(init, td, sf) endif endif endif return b[1:np] # optimize the (seasonal) initializer conditional on fixed @theta if sf == 2 matrix b = init[1:m+1] if verbose == 3 set max_verbose full elif verbose == 2 set max_verbose on set bfgs_verbskip 5 endif catch scalar SSE = BFGSmin(&b, ets_crit_mult(theta | b, my, m)) err = $error if err printf "Got error %d from BFGSmin\n", err return err endif init[1:m+2] = ets_mult_init(theta | b, m) else # the initializer is concentrated out by OLS ets_crit_add(theta, my, &init, td, m) endif if verbose > 0 printf "estimated initializer, given θ:\n" ets_print_init(init, td, sf) endif return 0 scalar td = b.trend scalar sf = b.seasonal scalar fc = ets_fcol(b) string s = sprintf("%s_level", b.yname) b[s] = "virtual_series:S[,2]" if td s = sprintf("%s_trend", b.yname) b[s] = "virtual_series:S[,3]" endif if sf > 0 s = sprintf("%s_season", b.yname) b[s] = sprintf("virtual_series:S[,%d]", 3 + td) endif s = sprintf("%s_fcast", b.yname) b[s] = sprintf("virtual_series:S[,%d]", fc) s = sprintf("%s_irreg", b.yname) b[s] = sprintf("virtual_series:S[,1]-S[,%d]", fc) scalar T = rows(b.S) scalar t0 = T - b.ntest + 1 scalar fc = ets_fcol(b) matrix e = b.S[t0:T,1] - b.S[t0:T,fc] b.RMSE_test = sqrt(sumc(e.^2) / b.ntest) b.MAE_test = sumc(abs(e)) / b.ntest # Forecast standard errors for the additive models (SES, Holt and # additive Holt-Winters), via the equivalent ETS(A,*,*) model: see # https://otexts.com/fpp3/ets-forecasting.html, section 8.7 and # Table 8.8. The h-step variance is # σ² [1 + Σ_{j=1}^{h-1} c_j²], # where c_j = α + β_e j + γ_e I(j mod m = 0). The ETS smoothing # parameters are related to ours by β_e = αβ and γ_e = (1-α)γ (see # fpp3 sections 8.3 and 8.5). scalar td = b.trend scalar m = b.seasonal > 0 ? rows(b.init) - 2 : 0 scalar α = b.theta[1] scalar β = td ? α * b.theta[2] : 0 scalar γ = m > 0 ? (1-α) * b.theta[3] : 0 scalar s2 = b.RMSE_train^2 scalar hmax = b.limits[3] - b.limits[2] matrix j = seq(1, hmax-1)' matrix cj = α + β * j if m > 0 cj += γ * ((j - m * floor(j/m)) .= 0) endif matrix vh = ones(rows(b.S), 1) * NA scalar T = rows(vh) matrix cs = hmax > 1 ? 0 | cum(cj.^2) : 0 vh[T-hmax+1:T] = s2 * (1 + cs) strings ss = cnameget(b.S) b.S ~= sqrt(vh) cnameset(b.S, ss + "stderr") # Simulated paths of y over the forecast range, multiplicative case, # for which there's no closed form for the forecast variance. Each # path runs the recursion forward from the final training state, # feeding in y = forecast + error, with the error drawn as # RMSE_train times Student t(df), so the interval at h = 1 matches # the analytic additive one. Returns an N x H matrix. scalar m = nelem(b.init) - 2 scalar H = b.limits[3] - b.limits[2] scalar n = rows(b.S) - H scalar α = b.theta[1] scalar β = b.theta[2] scalar γ = b.theta[3] matrix E = b.RMSE_train * mrandgen(t, df, N, H) matrix ι = ones(N, 1) matrix a_ = ι * b.S[n,2] matrix b_ = ι * b.S[n,3] # S[,h] holds s_{n+h-m}; S[,h+m] receives s_{n+h} matrix S = (ι * b.S[n-m+1:n,4]') ~ zeros(N, H) matrix Y = zeros(N, H) matrix a y sm loop h=1..H sm = S[,h] y = (a_ + b_) .* sm + E[,h] a = α * y ./ sm + (1-α) * (a_ + b_) b_ = β * (a - a_) + (1-β) * b_ S[,h+m] = γ * y ./ a + (1-γ) * sm a_ = a Y[,h] = y endloop return Y # For a c percent interval from the empirical quantiles of the # simulated paths @Y, return a T x 2 matrix holding the center and # half-width of the band, NA outside the forecast range (the last # cols(Y) of T rows). The band need not be symmetric about the # forecast, hence the representation via its midpoint. scalar H = cols(Y) scalar p = (100 - c) / 200 matrix Q = quantile(Y, {p, 1-p}) matrix ret = mshape(NA, T, 2) ret[T-H+1:T,1] = (Q[1,] + Q[2,])' / 2 ret[T-H+1:T,2] = (Q[2,] - Q[1,])' / 2 return ret scalar err = 0 if nelem(pc) != 2 err = 1 elif min(pc) < 50 || max(pc) > 99 err = 2 elif pc[1] == pc[2] err = 3 endif errorif(err, "ci_levels: invalid argument") if exists(ci_levels) ets_check_ci_levels(ci_levels) else ci_levels = {80, 95} endif hmax = b.limits[3] - b.limits[2] c1 = min(ci_levels) c2 = max(ci_levels) bundles bb = array(2) # for plotting, show the wider band first f1 = critical(t, df, (100 - c2) / 200) f2 = critical(t, df, (100 - c1) / 200) bb[1] = _(center="forecast", width="stderr", factor=f1, style="fill", color="#BBD9F5") bb[2] = _(center="forecast", width="stderr", factor=f2, style="fill", color="#6AADEC") return bb bundle b b.trend = 0 b.seasonal = 0 b.theta = {0.5}' # placeholder, length checked only if supplied b.init = {1.0} # ditto b.ntest = 0 b.verbose = 1 b.bounds = _(trend={0,1}, seasonal={0,2}, verbose={0,3}, ntest={0,NA}) return b # return non-zero if dataset is not suitable: we require # sufficient observations (seasonal models need more, and a # periodicity > 1, which is checked at run time) return $nobs < 10 bundle b b.alpha = _(depends="parms=2") b.beta = _(depends="parms=2") b.gamma = _(depends="parms=2") b.model = _(default="$pd>1 && $pd!=52 ? 4 : 2", maximum="$pd>1 && $pd<52 ? 4 : 2") return b include ETS.gfn # simple exponential smoothing: no trend, no seasonality open alg_exports.gdt --frompkg=ETS dataset addobs 5 bundle sb = ETS(Exports) ETS_plot(&sb, 1) # Holt: linear trend open auspop.gdt --frompkg=ETS dataset addobs 10 bundle hb = ETS(pop, _(trend=1)) ETS_plot(&hb, 1) ETS_plot(&hb, 2) # Holt-Winters, with 12 out-of-sample forecasts open aus_holidays.gdt --frompkg=ETS dataset addobs 12 bundle wb = ETS(Trips, _(seasonal=1)) ETS_plot(&wb, 1) ETS_plot(&wb, 2)