1 /*  R : A Computer Language for Statistical Data Analysis
2  *
3  *  Copyright (C) 2003-7  The R Core Team
4  *
5  *  This program is free software; you can redistribute it and/or modify
6  *  it under the terms of the GNU General Public License as published by
7  *  the Free Software Foundation; either version 2 of the License, or
8  *  (at your option) any later version.
9  *
10  *  This program is distributed in the hope that it will be useful,
11  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
12  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
13  *  GNU General Public License for more details.
14  *
15  *  You should have received a copy of the GNU General Public License
16  *  along with this program; if not, a copy is available at
17  *  https://www.R-project.org/Licenses/.
18  */
19 
20 /* Originally contributed by David Meyer */
21 
22 #include <stdlib.h>
23 #include <string.h>  // memcpy
24 
25 #include <R.h>
26 #include "ts.h"
27 
HoltWinters(double * x,int * xl,double * alpha,double * beta,double * gamma,int * start_time,int * seasonal,int * period,int * dotrend,int * doseasonal,double * a,double * b,double * s,double * SSE,double * level,double * trend,double * season)28 void HoltWinters (double *x,
29 		  int    *xl,
30 		  double *alpha,
31 		  double *beta,
32 		  double *gamma,
33 		  int    *start_time,
34 		  int    *seasonal,
35 		  int    *period,
36 		  int    *dotrend,
37 		  int    *doseasonal,
38 
39 		  double *a,
40 		  double *b,
41 		  double *s,
42 
43 		  /* return values */
44 		  double *SSE,
45 		  double *level,
46 		  double *trend,
47 		  double *season
48     )
49 
50 {
51     double res = 0, xhat = 0, stmp = 0;
52     int i, i0, s0;
53 
54     /* copy start values to the beginning of the vectors */
55     level[0] = *a;
56     if (*dotrend == 1) trend[0] = *b;
57     if (*doseasonal == 1) memcpy(season, s, *period * sizeof(double));
58 
59     for (i = *start_time - 1; i < *xl; i++) {
60 	/* indices for period i */
61 	i0 = i - *start_time + 2;
62 	s0 = i0 + *period - 1;
63 
64 	/* forecast *for* period i */
65 	xhat = level[i0 - 1] + (*dotrend == 1 ? trend[i0 - 1] : 0);
66 	stmp = *doseasonal == 1 ? season[s0 - *period] : (*seasonal != 1);
67 	if (*seasonal == 1)
68 	    xhat += stmp;
69 	else
70 	    xhat *= stmp;
71 
72 	/* Sum of Squared Errors */
73 	res   = x[i] - xhat;
74 	*SSE += res * res;
75 
76 	/* estimate of level *in* period i */
77 	if (*seasonal == 1)
78 	    level[i0] = *alpha       * (x[i] - stmp)
79 		      + (1 - *alpha) * (level[i0 - 1] + trend[i0 - 1]);
80 	else
81 	    level[i0] = *alpha       * (x[i] / stmp)
82 		      + (1 - *alpha) * (level[i0 - 1] + trend[i0 - 1]);
83 
84 	/* estimate of trend *in* period i */
85 	if (*dotrend == 1)
86 	    trend[i0] = *beta        * (level[i0] - level[i0 - 1])
87 		      + (1 - *beta)  * trend[i0 - 1];
88 
89 	/* estimate of seasonal component *in* period i */
90 	if (*doseasonal == 1) {
91 	    if (*seasonal == 1)
92 		season[s0] = *gamma       * (x[i] - level[i0])
93 			   + (1 - *gamma) * stmp;
94 	    else
95 		season[s0] = *gamma       * (x[i] / level[i0])
96 			   + (1 - *gamma) * stmp;
97 	}
98     }
99 }
100