timeseries.cpp
Go to the documentation of this file.
00001 /***************************************************************************
00002   file : $URL: https://frepple.svn.sourceforge.net/svnroot/frepple/trunk/modules/forecast/forecastsolver.cpp $
00003   version : $LastChangedRevision: 1656 $  $LastChangedBy: jdetaeye $
00004   date : $LastChangedDate: 2012-03-27 19:05:34 +0200 (Tue, 27 Mar 2012) $
00005  ***************************************************************************/
00006 
00007 /***************************************************************************
00008  *                                                                         *
00009  * Copyright (C) 2007-2012 by Johan De Taeye, frePPLe bvba                 *
00010  *                                                                         *
00011  * This library is free software; you can redistribute it and/or modify it *
00012  * under the terms of the GNU Lesser General Public License as published   *
00013  * by the Free Software Foundation; either version 2.1 of the License, or  *
00014  * (at your option) any later version.                                     *
00015  *                                                                         *
00016  * This library is distributed in the hope that it will be useful,         *
00017  * but WITHOUT ANY WARRANTY; without even the implied warranty of          *
00018  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser *
00019  * General Public License for more details.                                *
00020  *                                                                         *
00021  * You should have received a copy of the GNU Lesser General Public        *
00022  * License along with this library; if not, write to the Free Software     *
00023  * Foundation Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301,*
00024  * USA                                                                     *
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   // Validate the input
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   // Strip zero demand buckets at the start.
00047   // Eg when demand only starts in the middle of horizon, we only want to
00048   // use the second part of the horizon for forecasting. The zeros before the
00049   // demand start would distort the results.
00050   while (historycount >= 1 && *history == 0.0)
00051   {
00052     ++history;
00053     --historycount;
00054   }
00055 
00056   // We create the forecasting objects in stack memory for best performance.
00057   MovingAverage moving_avg;
00058   Croston croston;
00059   SingleExponential single_exp;
00060   DoubleExponential double_exp;
00061   Seasonal seasonal;
00062   int numberOfMethods = 4;
00063   ForecastMethod* methods[4];
00064 
00065   // Rules to determine which forecast methods can be applied
00066   methods[0] = &moving_avg;
00067   if (historycount < getForecastSkip() + 5)
00068     // Too little history: only use moving average
00069     numberOfMethods = 1;
00070   else
00071   {
00072     unsigned int zero = 0;
00073     for (unsigned long i = 0; i < historycount; ++i)
00074       if (!history[i]) ++zero;
00075     if (zero > Croston::getMinIntermittence() * historycount)
00076     {
00077       // If there are too many zeros: use croston or moving average.
00078       numberOfMethods = 2;
00079       methods[1] = &croston;
00080     }
00081     else
00082     {
00083       // The most common case: enough values and not intermittent
00084       methods[1] = &single_exp;
00085       methods[2] = &double_exp;
00086       methods[3] = &seasonal;
00087     }
00088   }
00089 
00090   // Initialize a vector with the smape weights
00091   double *weight = new double[historycount+1];
00092   weight[historycount] = 1.0;
00093   for (int i=historycount-1; i>=0; --i)
00094     weight[i] = weight[i+1] * Forecast::getForecastSmapeAlfa();
00095 
00096   // Evaluate each forecast method
00097   double best_error = DBL_MAX;
00098   int best_method = -1;
00099   double error;
00100   try
00101   {
00102     for (int i=0; i<numberOfMethods; ++i)
00103     {
00104       error = methods[i]->generateForecast(this, history, historycount, weight, debug);
00105       if (error<best_error)
00106       {
00107         best_error = error;
00108         best_method = i;
00109       }
00110     }
00111   }
00112   catch (...)
00113   {
00114     delete weight;
00115     throw;
00116   }
00117   delete weight;
00118 
00119   // Apply the most appropriate forecasting method
00120   if (best_method >= 0)
00121   {
00122     if (debug)
00123       logger << getName() << ": chosen method: " << methods[best_method]->getName() << endl;
00124     methods[best_method]->applyForecast(this, buckets, bucketcount, debug);
00125   }
00126 }
00127 
00128 
00129 //
00130 // MOVING AVERAGE FORECAST
00131 //
00132 
00133 
00134 unsigned int Forecast::MovingAverage::defaultbuckets = 5;
00135 
00136 
00137 double Forecast::MovingAverage::generateForecast
00138 (Forecast* fcst, const double history[], unsigned int count, const double weight[], bool debug)
00139 {
00140   double error_smape = 0.0;
00141 
00142   // Calculate the forecast and forecast error.
00143   for (unsigned int i = 1; i <= count; ++i)
00144   {
00145     double sum = 0.0;
00146     if (i > buckets)
00147     {
00148       for (unsigned int j = 0; j < buckets; ++j)
00149         sum += history[i-j-1];
00150       avg = sum / buckets;
00151     }
00152     else
00153     {
00154       // For the first few values
00155       for (unsigned int j = 0; j < i; ++j)
00156         sum += history[i-j-1];
00157       avg = sum / i;
00158     }
00159     if (i >= fcst->getForecastSkip() && i < count && avg + history[i] > ROUNDING_ERROR)
00160       error_smape += fabs(avg - history[i]) / (avg + history[i]) * weight[i];
00161   }
00162 
00163   // Echo the result
00164   if (debug)
00165     logger << (fcst ? fcst->getName() : "") << ": moving average : "
00166         << "smape " << error_smape
00167         << ", forecast " << avg << endl;
00168   return error_smape;
00169 }
00170 
00171 
00172 void Forecast::MovingAverage::applyForecast
00173 (Forecast* forecast, const Date buckets[], unsigned int bucketcount, bool debug)
00174 {
00175   // Loop over all buckets and set the forecast to a constant value
00176   if (avg < 0) return;
00177   for (unsigned int i = 1; i < bucketcount; ++i)
00178     forecast->setTotalQuantity(
00179       DateRange(buckets[i-1], buckets[i]),
00180       avg
00181     );
00182 }
00183 
00184 
00185 //
00186 // SINGLE EXPONENTIAL FORECAST
00187 //
00188 
00189 
00190 double Forecast::SingleExponential::initial_alfa = 0.2;
00191 double Forecast::SingleExponential::min_alfa = 0.03;
00192 double Forecast::SingleExponential::max_alfa = 1.0;
00193 
00194 
00195 double Forecast::SingleExponential::generateForecast
00196 (Forecast* fcst, const double history[], unsigned int count, const double weight[], bool debug)
00197 {
00198   // Verify whether this is a valid forecast method.
00199   //   - We need at least 5 buckets after the warmup period.
00200   if (count < fcst->getForecastSkip() + 5)
00201     return DBL_MAX;
00202 
00203   unsigned int iteration = 1;
00204   bool upperboundarytested = false;
00205   bool lowerboundarytested = false;
00206   double error = 0.0, error_smape = 0.0, best_smape = 0.0, delta, df_dalfa_i, sum_11, sum_12;
00207   double best_error = DBL_MAX, best_alfa = initial_alfa, best_f_i = 0.0;
00208   for (; iteration <= Forecast::getForecastIterations(); ++iteration)
00209   {
00210     // Initialize variables
00211     df_dalfa_i = sum_11 = sum_12 = error_smape = error = 0.0;
00212 
00213     // Initialize the iteration with the average of the first 3 values.
00214     f_i = (history[0] + history[1] + history[2]) / 3;
00215 
00216     // Calculate the forecast and forecast error.
00217     // We also compute the sums required for the Marquardt optimization.
00218     for (unsigned long i = 1; i <= count; ++i)
00219     {
00220       df_dalfa_i = history[i-1] - f_i + (1 - alfa) * df_dalfa_i;
00221       f_i = history[i-1] * alfa + (1 - alfa) * f_i;
00222       if (i == count) break;
00223       sum_12 += df_dalfa_i * (history[i] - f_i) * weight[i];
00224       sum_11 += df_dalfa_i * df_dalfa_i * weight[i];
00225       if (i >= fcst->getForecastSkip())
00226       {
00227         error += (f_i - history[i]) * (f_i - history[i]) * weight[i];
00228         if (f_i + history[i] > ROUNDING_ERROR)
00229           error_smape += fabs(f_i - history[i]) / (f_i + history[i]) * weight[i];
00230       }
00231     }
00232 
00233     // Better than earlier iterations?
00234     if (error < best_error)
00235     {
00236       best_error = error;
00237       best_smape = error_smape;
00238       best_alfa = alfa;
00239       best_f_i = f_i;
00240     }
00241 
00242     // Add Levenberg - Marquardt damping factor
00243     if (fabs(sum_11 + error / iteration) > ROUNDING_ERROR)
00244       sum_11 += error / iteration;
00245 
00246     // Calculate a delta for the alfa parameter
00247     if (fabs(sum_11) < ROUNDING_ERROR) break;
00248     delta = sum_12 / sum_11;
00249 
00250     // Stop when we are close enough and have tried hard enough
00251     if (fabs(delta) < ACCURACY && iteration > 3) break;
00252 
00253     // New alfa
00254     alfa += delta;
00255 
00256     // Stop when we repeatedly bounce against the limits.
00257     // Testing a limits once is allowed.
00258     if (alfa > max_alfa)
00259     {
00260       alfa = max_alfa;
00261       if (upperboundarytested) break;
00262       upperboundarytested = true;
00263     }
00264     else if (alfa < min_alfa)
00265     {
00266       alfa = min_alfa;
00267       if (lowerboundarytested) break;
00268       lowerboundarytested = true;
00269     }
00270   }
00271 
00272   // Keep the best result
00273   f_i = best_f_i;
00274 
00275   // Echo the result
00276   if (debug)
00277     logger << (fcst ? fcst->getName() : "") << ": single exponential : "
00278         << "alfa " << best_alfa
00279         << ", smape " << best_smape
00280         << ", " << iteration << " iterations"
00281         << ", forecast " << f_i << endl;
00282   return best_smape;
00283 }
00284 
00285 
00286 void Forecast::SingleExponential::applyForecast
00287 (Forecast* forecast, const Date buckets[], unsigned int bucketcount, bool debug)
00288 {
00289   // Loop over all buckets and set the forecast to a constant value
00290   if (f_i < 0) return;
00291   for (unsigned int i = 1; i < bucketcount; ++i)
00292     forecast->setTotalQuantity(
00293       DateRange(buckets[i-1], buckets[i]),
00294       f_i
00295     );
00296 }
00297 
00298 
00299 //
00300 // DOUBLE EXPONENTIAL FORECAST
00301 //
00302 
00303 double Forecast::DoubleExponential::initial_alfa = 0.2;
00304 double Forecast::DoubleExponential::min_alfa = 0.02;
00305 double Forecast::DoubleExponential::max_alfa = 1.0;
00306 double Forecast::DoubleExponential::initial_gamma = 0.2;
00307 double Forecast::DoubleExponential::min_gamma = 0.05;
00308 double Forecast::DoubleExponential::max_gamma = 1.0;
00309 double Forecast::DoubleExponential::dampenTrend = 0.8;
00310 
00311 
00312 double Forecast::DoubleExponential::generateForecast
00313 (Forecast* fcst, const double history[], unsigned int count, const double weight[], bool debug)
00314 {
00315   // Verify whether this is a valid forecast method.
00316   //   - We need at least 5 buckets after the warmup period.
00317   if (count < fcst->getForecastSkip() + 5)
00318     return DBL_MAX;
00319 
00320   // Define variables
00321   double error = 0.0, error_smape = 0.0, delta_alfa, delta_gamma, determinant;
00322   double constant_i_prev, trend_i_prev, d_constant_d_gamma_prev,
00323          d_constant_d_alfa_prev, d_constant_d_alfa, d_constant_d_gamma,
00324          d_trend_d_alfa, d_trend_d_gamma, d_forecast_d_alfa, d_forecast_d_gamma,
00325          sum11, sum12, sum22, sum13, sum23;
00326   double best_error = DBL_MAX, best_smape = 0, best_alfa = initial_alfa,
00327          best_gamma = initial_gamma, best_constant_i = 0.0, best_trend_i = 0.0;
00328 
00329   // Iterations
00330   unsigned int iteration = 1, boundarytested = 0;
00331   for (; iteration <= Forecast::getForecastIterations(); ++iteration)
00332   {
00333     // Initialize variables
00334     error = error_smape = sum11 = sum12 = sum22 = sum13 = sum23 = 0.0;
00335     d_constant_d_alfa = d_constant_d_gamma = d_trend_d_alfa = d_trend_d_gamma = 0.0;
00336     d_forecast_d_alfa = d_forecast_d_gamma = 0.0;
00337 
00338     // Initialize the iteration
00339     constant_i = (history[0] + history[1] + history[2]) / 3;
00340     trend_i = (history[3] - history[0]) / 3;
00341 
00342     // Calculate the forecast and forecast error.
00343     // We also compute the sums required for the Marquardt optimization.
00344     for (unsigned long i = 1; i <= count; ++i)
00345     {
00346       constant_i_prev = constant_i;
00347       trend_i_prev = trend_i;
00348       constant_i = history[i-1] * alfa + (1 - alfa) * (constant_i_prev + trend_i_prev);
00349       trend_i = gamma * (constant_i - constant_i_prev) + (1 - gamma) * trend_i_prev;
00350       if (i == count) break;
00351       d_constant_d_gamma_prev = d_constant_d_gamma;
00352       d_constant_d_alfa_prev = d_constant_d_alfa;
00353       d_constant_d_alfa = history[i-1] - constant_i_prev - trend_i_prev
00354           + (1 - alfa) * d_forecast_d_alfa;
00355       d_constant_d_gamma = (1 - alfa) * d_forecast_d_gamma;
00356       d_trend_d_alfa = gamma * (d_constant_d_alfa - d_constant_d_alfa_prev)
00357           + (1 - gamma) * d_trend_d_alfa;
00358       d_trend_d_gamma = constant_i - constant_i_prev - trend_i_prev
00359           + gamma * (d_constant_d_gamma - d_constant_d_gamma_prev)
00360           + (1 - gamma) * d_trend_d_gamma;
00361       d_forecast_d_alfa = d_constant_d_alfa + d_trend_d_alfa;
00362       d_forecast_d_gamma = d_constant_d_gamma + d_trend_d_gamma;
00363       sum11 += weight[i] * d_forecast_d_alfa * d_forecast_d_alfa;
00364       sum12 += weight[i] * d_forecast_d_alfa * d_forecast_d_gamma;
00365       sum22 += weight[i] * d_forecast_d_gamma * d_forecast_d_gamma;
00366       sum13 += weight[i] * d_forecast_d_alfa * (history[i] - constant_i - trend_i);
00367       sum23 += weight[i] * d_forecast_d_gamma * (history[i] - constant_i - trend_i);
00368       if (i >= fcst->getForecastSkip()) // Don't measure during the warmup period
00369       {
00370         error += (constant_i + trend_i - history[i]) * (constant_i + trend_i - history[i]) * weight[i];
00371         if (constant_i + trend_i + history[i] > ROUNDING_ERROR)
00372           error_smape += fabs(constant_i + trend_i - history[i]) / (constant_i + trend_i + history[i]) * weight[i];
00373       }
00374     }
00375 
00376     // Better than earlier iterations?
00377     if (error < best_error)
00378     {
00379       best_error = error;
00380       best_smape = error_smape;
00381       best_alfa = alfa;
00382       best_gamma = gamma;
00383       best_constant_i = constant_i;
00384       best_trend_i = trend_i;
00385     }
00386 
00387     // Add Levenberg - Marquardt damping factor
00388     sum11 += error / iteration;
00389     sum22 += error / iteration;
00390 
00391     // Calculate a delta for the alfa and gamma parameters
00392     determinant = sum11 * sum22 - sum12 * sum12;
00393     if (fabs(determinant) < ROUNDING_ERROR)
00394     {
00395       // Almost singular matrix. Try without the damping factor.
00396       sum11 -= error / iteration;
00397       sum22 -= error / iteration;
00398       determinant = sum11 * sum22 - sum12 * sum12;
00399       if (fabs(determinant) < ROUNDING_ERROR)
00400         // Still singular - stop iterations here
00401         break;
00402     }
00403     delta_alfa = (sum13 * sum22 - sum23 * sum12) / determinant;
00404     delta_gamma = (sum23 * sum11 - sum13 * sum12) / determinant;
00405 
00406     // Stop when we are close enough and have tried hard enough
00407     if (fabs(delta_alfa) + fabs(delta_gamma) < 2 * ACCURACY && iteration > 3)
00408       break;
00409 
00410     // New values for the next iteration
00411     alfa += delta_alfa;
00412     gamma += delta_gamma;
00413 
00414     // Limit the parameters in their allowed range.
00415     if (alfa > max_alfa)
00416       alfa = max_alfa;
00417     else if (alfa < min_alfa)
00418       alfa = min_alfa;
00419     if (gamma > max_gamma)
00420       gamma = max_gamma;
00421     else if (gamma < min_gamma)
00422       gamma = min_gamma;
00423 
00424     // Verify repeated running with both parameters at the boundary
00425     if ((gamma == min_gamma || gamma == max_gamma)
00426         && (alfa == min_alfa || alfa == max_alfa))
00427     {
00428       if (boundarytested++ > 5) break;
00429     }
00430   }
00431 
00432   // Keep the best result
00433   constant_i = best_constant_i;
00434   trend_i = best_trend_i;
00435 
00436   // Echo the result
00437   if (debug)
00438     logger << (fcst ? fcst->getName() : "") << ": double exponential : "
00439         << "alfa " << best_alfa
00440         << ", gamma " << best_gamma
00441         << ", smape " << best_smape
00442         << ", " << iteration << " iterations"
00443         << ", constant " << constant_i
00444         << ", trend " << trend_i
00445         << ", forecast " << (trend_i + constant_i) << endl;
00446   return best_smape;
00447 }
00448 
00449 
00450 void Forecast::DoubleExponential::applyForecast
00451 (Forecast* forecast, const Date buckets[], unsigned int bucketcount, bool debug)
00452 {
00453   // Loop over all buckets and set the forecast to a linearly changing value
00454   for (unsigned int i = 1; i < bucketcount; ++i)
00455   {
00456     constant_i += trend_i;
00457     trend_i *= dampenTrend; // Reduce slope in the future
00458     if (constant_i > 0)
00459       forecast->setTotalQuantity(
00460         DateRange(buckets[i-1], buckets[i]),
00461         constant_i
00462       );
00463   }
00464 }
00465 
00466 
00467 //
00468 // SEASONAL FORECAST
00469 //
00470 
00471 unsigned int Forecast::Seasonal::min_period = 2;
00472 unsigned int Forecast::Seasonal::max_period = 14;
00473 double Forecast::Seasonal::dampenTrend = 0.8;
00474 double Forecast::Seasonal::initial_alfa = 0.2;
00475 double Forecast::Seasonal::min_alfa = 0.02;
00476 double Forecast::Seasonal::max_alfa = 1.0;
00477 double Forecast::Seasonal::initial_beta = 0.2;
00478 double Forecast::Seasonal::min_beta = 0.2;
00479 double Forecast::Seasonal::max_beta = 1.0;
00480 double Forecast::Seasonal::initial_gamma = 0.3;
00481 double Forecast::Seasonal::min_gamma = 0.05;
00482 double Forecast::Seasonal::max_gamma = 1.0;
00483 
00484 
00485 void Forecast::Seasonal::detectCycle(const double history[], unsigned int count)
00486 {
00487   // We need at least 2 cycles
00488   if (count < min_period*2) return;
00489 
00490   // Compute the average value
00491   double average = 0.0;
00492   for (unsigned int i = 0; i < count; ++i)
00493     average += history[i];
00494   average /= count;
00495 
00496   // Compute variance
00497   double variance = 0.0;
00498   for (unsigned int i = 0; i < count; ++i)
00499     variance += (history[i]-average)*(history[i]-average);
00500   variance /= count;
00501 
00502   // Compute autocorrelation for different periods
00503   double prev = 10.0;
00504   double prevprev = 10.0;
00505   for (unsigned short p = min_period; p <= max_period && p < count/2; ++p)
00506   {
00507     // Compute correlation
00508     double correlation = 0.0;
00509     for (unsigned int i = p; i < count; ++i)
00510       correlation += (history[i]-average)*(history[i-p]-average);
00511     correlation /= count - p;
00512     correlation /= variance;
00513     // Detect cycles if we find a period with a big autocorrelation which
00514     // is significantly larger than the adjacent periods.
00515     if ((p > min_period + 1 && prev > prevprev*1.1) && prev > correlation*1.1 && prev > 0.3)
00516     {
00517       period = p - 1;
00518       return;
00519     }
00520     prevprev = prev;
00521     prev = correlation;
00522   }
00523 }
00524 
00525 
00526 double Forecast::Seasonal::generateForecast
00527 (Forecast* fcst, const double history[], unsigned int count, const double weight[], bool debug)
00528 {
00529   // Check for seasonal cycles
00530   detectCycle(history, count);
00531 
00532   // Return if no seasonality is found
00533   if (!period) return DBL_MAX;
00534 
00535   // Define variables
00536   double error = 0.0, error_smape = 0.0, delta_alfa, delta_beta, delta_gamma;
00537   double sum11, sum12, sum13, sum14, sum22, sum23, sum24, sum33, sum34;
00538   double best_error = DBL_MAX, best_smape = 0, best_alfa = initial_alfa,
00539          best_beta = initial_beta, best_gamma = initial_gamma;
00540   S_i = new double[period];
00541 
00542   // Iterations
00543   double L_i_prev;
00544   unsigned int iteration = 1, boundarytested = 0;
00545   for (; iteration <= Forecast::getForecastIterations(); ++iteration)
00546   {
00547     // Initialize variables
00548     error = error_smape = sum11 = sum12 = sum13 = sum14 = 0.0;
00549     sum22 = sum23 = sum24 = sum33 = sum34 = 0.0;
00550 
00551     // Initialize the iteration
00552     // L_i = average over first cycle
00553     // T_i = average delta measured in second cycle
00554     // S_i[index] = actual divided by average in first cycle
00555     L_i = 0.0;
00556     for (unsigned long i = 0; i < period; ++i)
00557       L_i += history[i];
00558     L_i /= period;
00559     T_i = 0.0;
00560     for (unsigned long i = 0; i < period; ++i)
00561     {
00562       T_i += history[i+period] - history[i];
00563       S_i[i] = history[i] / L_i;
00564     }
00565     T_i /= period * period;
00566 
00567     // Calculate the forecast and forecast error.
00568     // We also compute the sums required for the Marquardt optimization.
00569     cycleindex = 0;
00570     for (unsigned long i = period; i <= count; ++i)
00571     {
00572       L_i_prev = L_i;
00573       if (S_i[cycleindex] > ROUNDING_ERROR)
00574         L_i = alfa * history[i-1] / S_i[cycleindex] + (1 - alfa) * (L_i + T_i);
00575       else
00576         L_i = (1 - alfa) * (L_i + T_i);
00577       T_i = beta * (L_i - L_i_prev) + (1 - beta) * T_i;
00578       S_i[cycleindex] = gamma * history[i-1] / L_i + (1 - gamma) * S_i[cycleindex];
00579       if (i == count) break;
00580       if (i >= fcst->getForecastSkip()) // Don't measure during the warmup period
00581       {
00582         double fcst = (L_i + T_i) * S_i[cycleindex];
00583         error += (fcst - history[i]) * (fcst - history[i]) * weight[i];
00584         if (fcst + history[i] > ROUNDING_ERROR)
00585           error_smape += fabs(fcst - history[i]) / (fcst + history[i]) * weight[i];
00586       }
00587       if (++cycleindex >= period) cycleindex = 0;
00588     }
00589 
00590     // Better than earlier iterations?
00591     best_smape = error_smape;
00592     if (error < best_error)
00593     {
00594       best_error = error;
00595       best_smape = error_smape;
00596       best_alfa = alfa;
00597       best_beta = beta;
00598       best_gamma = gamma;
00599     }
00600     break; // @todo no iterations yet to tune the seasonal parameters
00601 
00602     // Add Levenberg - Marquardt damping factor
00603     sum11 += error / iteration;
00604     sum22 += error / iteration;
00605     sum33 += error / iteration;
00606 
00607     // Calculate a delta for the alfa, beta and gamma parameters.
00608     // We're using Cramer's rule to solve a set of linear equations.
00609     double det = determinant(sum11, sum12, sum13,
00610         sum12, sum22, sum23,
00611         sum13, sum23, sum33);
00612     if (fabs(det) < ROUNDING_ERROR)
00613     {
00614       // Almost singular matrix. Try without the damping factor.
00615       sum11 -= error / iteration;
00616       sum22 -= error / iteration;
00617       sum33 -= error / iteration;
00618       det = determinant(sum11, sum12, sum13,
00619           sum12, sum22, sum23,
00620           sum13, sum23, sum33);
00621       if (fabs(det) < ROUNDING_ERROR)
00622         // Still singular - stop iterations here
00623         break;
00624     }
00625     delta_alfa = determinant(sum14, sum12, sum13,
00626         sum24, sum22, sum23,
00627         sum34, sum23, sum33) / det;
00628     delta_beta = determinant(sum11, sum14, sum13,
00629         sum12, sum24, sum23,
00630         sum13, sum34, sum33) / det;
00631     delta_gamma = determinant(sum11, sum12, sum14,
00632         sum12, sum22, sum24,
00633         sum13, sum23, sum34) / det;
00634 
00635     // Stop when we are close enough and have tried hard enough
00636     if (fabs(delta_alfa) + fabs(delta_beta) + fabs(delta_gamma) < 3 * ACCURACY
00637         && iteration > 3)
00638       break;
00639 
00640     // New values for the next iteration
00641     alfa += delta_alfa;
00642     alfa += delta_beta;
00643     gamma += delta_gamma;
00644 
00645     // Limit the parameters in their allowed range.
00646     if (alfa > max_alfa)
00647       alfa = max_alfa;
00648     else if (alfa < min_alfa)
00649       alfa = min_alfa;
00650     if (beta > max_beta)
00651       beta = max_beta;
00652     else if (beta < min_beta)
00653       beta = min_beta;
00654     if (gamma > max_gamma)
00655       gamma = max_gamma;
00656     else if (gamma < min_gamma)
00657       gamma = min_gamma;
00658 
00659     // Verify repeated running with any parameters at the boundary
00660     if (gamma == min_gamma || gamma == max_gamma ||
00661         beta == min_beta || beta == max_beta ||
00662         alfa == min_alfa || alfa == max_alfa)
00663     {
00664       if (boundarytested++ > 5) break;
00665     }
00666   }
00667 
00668   if (period > fcst->getForecastSkip())
00669     // Correction on the error: we counted less buckets. We now
00670     // proportionally increase the error to account for this and have a
00671     // value that can be compared with the other forecast methods.
00672     best_smape *= (count - fcst->getForecastSkip()) / (count - period);
00673 
00674   // Keep the best result
00675 
00676   // Echo the result
00677   if (debug)
00678     logger << (fcst ? fcst->getName() : "") << ": seasonal : "
00679         << "alfa " << best_alfa
00680         << ", beta " << best_beta
00681         << ", gamma " << best_gamma
00682         << ", smape " << best_smape
00683         << ", " << iteration << " iterations"
00684         << ", period " << period
00685         << ", constant " << L_i
00686         << ", trend " << T_i
00687         << ", forecast " << (L_i + T_i) * S_i[count % period]
00688         << endl;
00689   return best_smape;
00690 }
00691 
00692 
00693 void Forecast::Seasonal::applyForecast
00694 (Forecast* forecast, const Date buckets[], unsigned int bucketcount, bool debug)
00695 {
00696   // Loop over all buckets and set the forecast to a linearly changing value
00697   for (unsigned int i = 1; i < bucketcount; ++i)
00698   {
00699     L_i += T_i;
00700     T_i *= dampenTrend; // Reduce slope in the future
00701     if (L_i * S_i[cycleindex] > 0)
00702       forecast->setTotalQuantity(
00703         DateRange(buckets[i-1], buckets[i]),
00704         L_i * S_i[cycleindex]
00705       );
00706     if (++cycleindex >= period) cycleindex = 0;
00707   }
00708 }
00709 
00710 
00711 //
00712 // CROSTON'S FORECAST METHOD
00713 //
00714 
00715 
00716 double Forecast::Croston::initial_alfa = 0.1;
00717 double Forecast::Croston::min_alfa = 0.03;
00718 double Forecast::Croston::max_alfa = 1.0;
00719 double Forecast::Croston::min_intermittence = 0.33;
00720 
00721 
00722 double Forecast::Croston::generateForecast
00723 (Forecast* fcst, const double history[], unsigned int count, const double weight[], bool debug)
00724 {
00725   unsigned int iteration = 1;
00726   bool upperboundarytested = false;
00727   bool lowerboundarytested = false;
00728   double error = 0.0, error_smape = 0.0, best_smape = 0.0, delta;
00729   double q_i, p_i, d_p_i, d_q_i, d_f_i, sum1, sum2;
00730   double best_error = DBL_MAX, best_alfa = initial_alfa, best_f_i = 0.0;
00731   unsigned int between_demands = 1;
00732   for (; iteration <= Forecast::getForecastIterations(); ++iteration)
00733   {
00734     // Initialize variables
00735     error_smape = error = d_p_i = d_q_i = d_f_i = sum1 = sum2 = 0.0;
00736 
00737     // Initialize the iteration.
00738     q_i = f_i = history[0];
00739     p_i = 0;
00740 
00741     // Calculate the forecast and forecast error.
00742     // We also compute the sums required for the Marquardt optimization.
00743     for (unsigned long i = 1; i <= count; ++i)
00744     {
00745       if (history[i-1])
00746       {
00747         // Non-zero bucket
00748         d_p_i = between_demands - p_i + (1 - alfa) * d_p_i;
00749         d_q_i = history[i-1] - q_i + (1 - alfa) * d_q_i;
00750         q_i = alfa * history[i-1] + (1 - alfa) * q_i;
00751         p_i = alfa * between_demands + (1 - alfa) * q_i;
00752         f_i = q_i / p_i;
00753         d_f_i = (d_q_i - d_p_i * q_i / p_i) / p_i;
00754         between_demands = 1;
00755       }
00756       else
00757         ++between_demands;
00758       if (i == count) break;
00759       sum1 += weight[i] * d_f_i * (history[i] - f_i);
00760       sum2 += weight[i] * d_f_i * d_f_i;
00761       if (i >= fcst->getForecastSkip() && p_i > 0)
00762       {
00763         error += (f_i - history[i]) * (f_i - history[i]) * weight[i];
00764         if (f_i + history[i] > ROUNDING_ERROR)
00765           error_smape += fabs(f_i - history[i]) / (f_i + history[i]) * weight[i];
00766       }
00767     }
00768 
00769     // Better than earlier iterations?
00770     if (error < best_error)
00771     {
00772       best_error = error;
00773       best_smape = error_smape;
00774       best_alfa = alfa;
00775       best_f_i = f_i;
00776     }
00777 
00778     // Add Levenberg - Marquardt damping factor
00779     if (fabs(sum2 + error / iteration) > ROUNDING_ERROR)
00780       sum2 += error / iteration;
00781 
00782     // Calculate a delta for the alfa parameter
00783     if (fabs(sum2) < ROUNDING_ERROR) break;
00784     delta = sum1 / sum2;
00785 
00786     // Stop when we are close enough and have tried hard enough
00787     if (fabs(delta) < ACCURACY && iteration > 3) break;
00788 
00789     // New alfa
00790     alfa += delta;
00791 
00792     // Stop when we repeatedly bounce against the limits.
00793     // Testing a limits once is allowed.
00794     if (alfa > max_alfa)
00795     {
00796       alfa = max_alfa;
00797       if (upperboundarytested) break;
00798       upperboundarytested = true;
00799     }
00800     else if (alfa < min_alfa)
00801     {
00802       alfa = min_alfa;
00803       if (lowerboundarytested) break;
00804       lowerboundarytested = true;
00805     }
00806   }
00807 
00808   // Keep the best result
00809   f_i = best_f_i;
00810   alfa = best_alfa;
00811 
00812   // Echo the result
00813   if (debug)
00814     logger << (fcst ? fcst->getName() : "") << ": croston : "
00815         << "alfa " << best_alfa
00816         << ", smape " << best_smape
00817         << ", " << iteration << " iterations"
00818         << ", forecast " << f_i << endl;
00819   return best_smape;
00820 }
00821 
00822 
00823 void Forecast::Croston::applyForecast
00824 (Forecast* forecast, const Date buckets[], unsigned int bucketcount, bool debug)
00825 {
00826   // Loop over all buckets and set the forecast to a constant value
00827   if (f_i < 0) return;
00828   for (unsigned int i = 1; i < bucketcount; ++i)
00829     forecast->setTotalQuantity(
00830       DateRange(buckets[i-1], buckets[i]),
00831       f_i
00832     );
00833 }
00834 
00835 }       // end namespace

Documentation generated for frePPLe by  doxygen