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