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

Documentation generated for frePPLe by  doxygen