00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028 #include "forecast.h"
00029
00030 namespace module_forecast
00031 {
00032
00033 #define ACCURACY 0.01
00034
00035 void Forecast::generateFutureValues(
00036 const double history[], unsigned int historycount,
00037 const Date buckets[], unsigned int bucketcount,
00038 bool debug)
00039 {
00040
00041 if (!history || !buckets)
00042 throw RuntimeException("Null argument to forecast function");
00043 if (bucketcount < 2)
00044 throw DataException("Need at least 2 forecast dates");
00045
00046
00047
00048 MovingAverage moving_avg;
00049 SingleExponential single_exp;
00050 DoubleExponential double_exp;
00051 int numberOfMethods = 3;
00052 ForecastMethod* methods[3];
00053 methods[0] = &single_exp;
00054 methods[1] = &double_exp;
00055 methods[2] = &moving_avg;
00056
00057
00058 double *madWeight = new double[historycount+1];
00059 madWeight[historycount] = madWeight[historycount-1] = 1.0;
00060 for (int i=historycount-2; i>=0; --i)
00061 madWeight[i] = madWeight[i+1] * Forecast::getForecastMadAlfa();
00062
00063
00064 double best_error = DBL_MAX;
00065 int best_method = -1;
00066 double error;
00067 try
00068 {
00069 for (int i=0; i<numberOfMethods; ++i)
00070 {
00071 error = methods[i]->generateForecast(this, history, historycount, madWeight, debug);
00072 if (error<best_error)
00073 {
00074 best_error = error;
00075 best_method = i;
00076 }
00077 }
00078 }
00079 catch (...)
00080 {
00081 delete madWeight;
00082 throw;
00083 }
00084 delete madWeight;
00085
00086
00087 if (best_method >= 0)
00088 {
00089 if (debug)
00090 logger << getName() << ": chosen method : " << methods[best_method]->getName() << endl;
00091 methods[best_method]->applyForecast(this, buckets, bucketcount, debug);
00092 }
00093 }
00094
00095
00096
00097
00098
00099
00100
00101 unsigned int Forecast::MovingAverage::defaultbuckets = 5;
00102
00103
00104 double Forecast::MovingAverage::generateForecast
00105 (Forecast* fcst, const double history[], unsigned int count, const double madWeight[], bool debug)
00106 {
00107 double error_mad = 0.0;
00108 double sum = 0.0;
00109
00110
00111 for (unsigned int i = 1; i <= count; ++i)
00112 {
00113 sum += history[i-1];
00114 if (i>buckets)
00115 {
00116 sum -= history[i-1-buckets];
00117 avg = sum / buckets;
00118 }
00119 else
00120 avg = sum / i;
00121 if (i >= fcst->getForecastSkip())
00122 error_mad += fabs(avg - history[i]) * madWeight[i];
00123 }
00124
00125
00126 if (debug)
00127 logger << (fcst ? fcst->getName() : "") << ": moving average : "
00128 << "mad " << error_mad
00129 << ", forecast " << avg << endl;
00130 return error_mad;
00131 }
00132
00133
00134 void Forecast::MovingAverage::applyForecast
00135 (Forecast* forecast, const Date buckets[], unsigned int bucketcount, bool debug)
00136 {
00137
00138 if (avg < 0) return;
00139 for (unsigned int i = 1; i < bucketcount; ++i)
00140 forecast->setTotalQuantity(
00141 DateRange(buckets[i-1], buckets[i]),
00142 avg
00143 );
00144 }
00145
00146
00147
00148
00149
00150
00151
00152 double Forecast::SingleExponential::initial_alfa = 0.2;
00153 double Forecast::SingleExponential::min_alfa = 0.03;
00154 double Forecast::SingleExponential::max_alfa = 1.0;
00155
00156
00157 double Forecast::SingleExponential::generateForecast
00158 (Forecast* fcst, const double history[], unsigned int count, const double madWeight[], bool debug)
00159 {
00160
00161
00162 if (count < fcst->getForecastSkip() + 5)
00163 return DBL_MAX;
00164
00165 unsigned int iteration = 1;
00166 bool upperboundarytested = false;
00167 bool lowerboundarytested = false;
00168 double error_mad = 0.0, delta, df_dalfa_i, sum_11, sum_12;
00169 double best_error_mad = DBL_MAX, best_alfa = initial_alfa, best_f_i = 0.0;
00170 for (; iteration <= Forecast::getForecastIterations(); ++iteration)
00171 {
00172
00173 f_i = history[0];
00174 df_dalfa_i = sum_11 = sum_12 = error_mad = 0.0;
00175
00176
00177
00178 for (unsigned long i = 1; i <= count; ++i)
00179 {
00180 df_dalfa_i = history[i-1] - f_i + (1 - alfa) * df_dalfa_i;
00181 f_i = history[i-1] * alfa + (1 - alfa) * f_i;
00182 if (i == count) break;
00183 sum_12 += df_dalfa_i * (history[i] - f_i) * madWeight[i];
00184 sum_11 += df_dalfa_i * df_dalfa_i * madWeight[i];
00185 if (i >= fcst->getForecastSkip())
00186 error_mad += fabs(f_i - history[i]) * madWeight[i];
00187 }
00188
00189
00190 if (error_mad < best_error_mad)
00191 {
00192 best_error_mad = error_mad;
00193 best_alfa = alfa;
00194 best_f_i = f_i;
00195 }
00196
00197
00198 if (fabs(sum_11 + error_mad / iteration) > ROUNDING_ERROR)
00199 sum_11 += error_mad / iteration;
00200
00201
00202 if (fabs(sum_11) < ROUNDING_ERROR) break;
00203 delta = sum_12 / sum_11;
00204
00205
00206 if (fabs(delta) < ACCURACY && iteration > 3) break;
00207
00208
00209 alfa += delta;
00210
00211
00212
00213 if (alfa > max_alfa)
00214 {
00215 alfa = max_alfa;
00216 if (upperboundarytested) break;
00217 upperboundarytested = true;
00218 }
00219 else if (alfa < min_alfa)
00220 {
00221 alfa = min_alfa;
00222 if (lowerboundarytested) break;
00223 lowerboundarytested = true;
00224 }
00225 }
00226
00227
00228 f_i = best_f_i;
00229
00230
00231 if (debug)
00232 logger << (fcst ? fcst->getName() : "") << ": single exponential : "
00233 << "alfa " << best_alfa
00234 << ", mad " << best_error_mad
00235 << ", " << iteration << " iterations"
00236 << ", forecast " << f_i << endl;
00237 return best_error_mad;
00238 }
00239
00240
00241 void Forecast::SingleExponential::applyForecast
00242 (Forecast* forecast, const Date buckets[], unsigned int bucketcount, bool debug)
00243 {
00244
00245 if (f_i < 0) return;
00246 for (unsigned int i = 1; i < bucketcount; ++i)
00247 forecast->setTotalQuantity(
00248 DateRange(buckets[i-1], buckets[i]),
00249 f_i
00250 );
00251 }
00252
00253
00254
00255
00256
00257
00258
00259 double Forecast::DoubleExponential::initial_alfa = 0.2;
00260 double Forecast::DoubleExponential::min_alfa = 0.02;
00261 double Forecast::DoubleExponential::max_alfa = 1.0;
00262 double Forecast::DoubleExponential::initial_gamma = 0.2;
00263 double Forecast::DoubleExponential::min_gamma = 0.05;
00264 double Forecast::DoubleExponential::max_gamma = 1.0;
00265
00266
00267 double Forecast::DoubleExponential::generateForecast
00268 (Forecast* fcst, const double history[], unsigned int count, const double madWeight[], bool debug)
00269 {
00270
00271
00272 if (count < fcst->getForecastSkip() + 5)
00273 return DBL_MAX;
00274
00275
00276 double error_mad = 0.0, delta_alfa, delta_gamma, determinant;
00277 double constant_i_prev, trend_i_prev, d_constant_d_gamma_prev,
00278 d_constant_d_alfa_prev, d_constant_d_alfa, d_constant_d_gamma,
00279 d_trend_d_alfa, d_trend_d_gamma, d_forecast_d_alfa, d_forecast_d_gamma,
00280 sum11, sum12, sum22, sum13, sum23;
00281 double best_error_mad = DBL_MAX, best_alfa = initial_alfa,
00282 best_gamma = initial_gamma, best_constant_i = 0.0, best_trend_i = 0.0;
00283
00284
00285 unsigned int iteration = 1, boundarytested = 0;
00286 for (; iteration <= Forecast::getForecastIterations(); ++iteration)
00287 {
00288
00289 constant_i = history[0];
00290 trend_i = history[1] - history[0];
00291 error_mad = sum11 = sum12 = sum22 = sum13 = sum23 = 0.0;
00292 d_constant_d_alfa = d_constant_d_gamma = d_trend_d_alfa = d_trend_d_gamma = 0.0;
00293 d_forecast_d_alfa = d_forecast_d_gamma = 0.0;
00294
00295
00296
00297 for (unsigned long i = 1; i <= count; ++i)
00298 {
00299 constant_i_prev = constant_i;
00300 trend_i_prev = trend_i;
00301 constant_i = history[i-1] * alfa + (1 - alfa) * (constant_i_prev + trend_i_prev);
00302 trend_i = gamma * (constant_i - constant_i_prev) + (1 - gamma) * trend_i_prev;
00303 if (i == count) break;
00304 d_constant_d_gamma_prev = d_constant_d_gamma;
00305 d_constant_d_alfa_prev = d_constant_d_alfa;
00306 d_constant_d_alfa = history[i-1] - constant_i_prev - trend_i_prev
00307 + (1 - alfa) * d_forecast_d_alfa;
00308 d_constant_d_gamma = (1 - alfa) * d_forecast_d_gamma;
00309 d_trend_d_alfa = gamma * (d_constant_d_alfa - d_constant_d_alfa_prev)
00310 + (1 - gamma) * d_trend_d_alfa;
00311 d_trend_d_gamma = constant_i - constant_i_prev - trend_i_prev
00312 + gamma * (d_constant_d_gamma - d_constant_d_gamma_prev)
00313 + (1 - gamma) * d_trend_d_gamma;
00314 d_forecast_d_alfa = d_constant_d_alfa + d_trend_d_alfa;
00315 d_forecast_d_gamma = d_constant_d_gamma + d_trend_d_gamma;
00316 sum11 += madWeight[i] * d_forecast_d_alfa * d_forecast_d_alfa;
00317 sum12 += madWeight[i] * d_forecast_d_alfa * d_forecast_d_gamma;
00318 sum22 += madWeight[i] * d_forecast_d_gamma * d_forecast_d_gamma;
00319 sum13 += madWeight[i] * d_forecast_d_alfa * (history[i] - constant_i - trend_i);
00320 sum23 += madWeight[i] * d_forecast_d_gamma * (history[i] - constant_i - trend_i);
00321 if (i >= fcst->getForecastSkip())
00322 error_mad += fabs(constant_i + trend_i - history[i]) * madWeight[i];
00323 }
00324
00325
00326 if (error_mad < best_error_mad)
00327 {
00328 best_error_mad = error_mad;
00329 best_alfa = alfa;
00330 best_gamma = gamma;
00331 best_constant_i = constant_i;
00332 best_trend_i = trend_i;
00333 }
00334
00335
00336 sum11 += error_mad / iteration;
00337 sum22 += error_mad / iteration;
00338
00339
00340 determinant = sum11 * sum22 - sum12 * sum12;
00341 if (fabs(determinant) < ROUNDING_ERROR)
00342 {
00343
00344 sum11 -= error_mad / iteration;
00345 sum22 -= error_mad / iteration;
00346 determinant = sum11 * sum22 - sum12 * sum12;
00347 }
00348 if (fabs(determinant) < ROUNDING_ERROR) break;
00349 delta_alfa = (sum13 * sum22 - sum23 * sum12) / determinant;
00350 delta_gamma = (sum23 * sum11 - sum13 * sum12) / determinant;
00351
00352
00353 if (fabs(delta_alfa) + fabs(delta_gamma) < ACCURACY && iteration > 3)
00354 break;
00355
00356
00357 alfa += delta_alfa;
00358 gamma += delta_gamma;
00359
00360
00361 if (alfa > max_alfa)
00362 alfa = max_alfa;
00363 else if (alfa < min_alfa)
00364 alfa = min_alfa;
00365 if (gamma > max_gamma)
00366 gamma = max_gamma;
00367 else if (gamma < min_gamma)
00368 gamma = min_gamma;
00369
00370
00371 if ((gamma == min_gamma || gamma == max_gamma)
00372 && (alfa == min_alfa || alfa == max_alfa))
00373 {
00374 if (boundarytested++ > 2) break;
00375 }
00376 }
00377
00378
00379 constant_i = best_constant_i;
00380 trend_i = best_trend_i;
00381
00382
00383 if (debug)
00384 logger << (fcst ? fcst->getName() : "") << ": double exponential : "
00385 << "alfa " << best_alfa
00386 << ", gamma " << best_gamma
00387 << ", mad " << best_error_mad
00388 << ", " << iteration << " iterations"
00389 << ", constant " << constant_i
00390 << ", trend " << trend_i
00391 << ", forecast " << (trend_i + constant_i) << endl;
00392 return best_error_mad;
00393 }
00394
00395
00396 void Forecast::DoubleExponential::applyForecast
00397 (Forecast* forecast, const Date buckets[], unsigned int bucketcount, bool debug)
00398 {
00399
00400 for (unsigned int i = 1; i < bucketcount; ++i)
00401 if (constant_i + i * trend_i > 0)
00402 forecast->setTotalQuantity(
00403 DateRange(buckets[i-1], buckets[i]),
00404 constant_i + i * trend_i
00405 );
00406 }
00407
00408
00409 }