I created a new Robust HoltWinters function (based on the stats :: Holt-Winters method) in R (for "Strong forecasting with exponential and hilly winters" Smoothing "by Sarah Gelper1, Roland Freed, Christoph Crox, September 26, 2008) Why? Well? ... why not! But I digress ...
Kernel of statistics: the Holt-Winters method is C code called C_HoltWinters , which I modified to be reliable (see below)
#include <stdlib.h> #include <string.h> // memcpy #include <math.h> #include <Rh> #include "ts.h" void HoltWinters ( double *x, /*as.double(x) */ double *x_adj, /*Adjust time series data, if need be Added*/ int *xl, /*lenx - Length of the current time series*/ double *alpha, /*as.double(max(min(alpha, 1), 0)), */ double *beta, /*as.double(max(min(beta,1), 0)), */ double *gamma, /*as.double(max(min(gamma, 1), 0)), */ double *llamda,/*as.double(max(min(llamda,1),0)), ADDED*/ int *start_time, /*as.integer(start.time), */ int *seasonal, /*as.integer(!+(seasonal == "multiplicative")), */ int *period, /* as.integer(f), */ int *dotrend, /* as.integer(!is.logical(beta) || beta), */ int *doseasonal, /* as.integer(!is.logical(gamma) || gamma), */ double *a, /*l.start - starting values for level*/ double *b, /*b.start - starting values for Trend*/ double *s, /*s.start - starting values for SEasonal*/ double *l, /*t.start - starting values for LLamda ADDED*/ double *k, /* Value for K ADDED*/ double *ck, /*value for ck ADDED*/ /* return values */ double *SSE, double *level, double *trend, double *season ) { double res = 0, xhat = 0, stmp = 0, theta = 1, RhoK = 0, phi = 0 ; int i, i0, s0; /*i is the current t, i0 is the current LESS starting period, and s0 = is the seasonal current LESS Starting period*/ /* copy start values to the beginning of the vectors */ level[0] = *a; if (*dotrend == 1) trend[0] = *b; if (*doseasonal == 1) memcpy(season, s, *period * sizeof(double)); for (i = *start_time - 1; i < *xl; i++) { /* indices for period i */ i0 = i - *start_time + 2; s0 = i0 + *period - 1; /* forecast *for* period i */ xhat = level[i0 - 1] + (*dotrend == 1 ? trend[i0 - 1] : 0); stmp = *doseasonal == 1 ? season[s0 - *period] : (*seasonal != 1); if (*seasonal == 1) xhat += stmp; else xhat *= stmp; /* Sum of Squared Errors */ res = x[i] - xhat; /*adjusting for robustness....Gahds*/ RhoK = (abs(res / theta) <= *k ? *ck * (1 - pow(1 - pow((res / (*k * theta)),2),3)): *ck); theta = sqrt(*llamda * RhoK * pow(theta,2) + (1 - *llamda) * pow(theta,2)); phi = (abs(res / theta) < *k ? res / theta : ((res / theta) / abs(res / theta) * (*k))); x_adj[i] = phi * theta + xhat; res = x_adj[i] - xhat; *SSE += res * res; /* estimate of level *in* period i */ if (*seasonal == 1) level[i0] = *alpha * (x_adj[i] - stmp) + (1 - *alpha) * (level[i0 - 1] + trend[i0 - 1]); else level[i0] = *alpha * (x_adj[i] / stmp) + (1 - *alpha) * (level[i0 - 1] + trend[i0 - 1]); /* estimate of trend *in* period i */ if (*dotrend == 1) trend[i0] = *beta * (level[i0] - level[i0 - 1]) + (1 - *beta) * trend[i0 - 1]; /* estimate of seasonal component *in* period i */ if (*doseasonal == 1) { if (*seasonal == 1) season[s0] = *gamma * (x_adj[i] - level[i0]) + (1 - *gamma) * stmp; else season[s0] = *gamma * (x_adj[i] / level[i0]) + (1 - *gamma) * stmp; } } }
So, I will compile it in a sigh of Windows using R (3.2.2):
R CMD SHLIB C_R_HoltWinters.c gcc -m64 -I"C:/PROGRA~1/R/R-32~1.2/include" -DNDEBUG -I"d:/RCompile/r-compiling/local/local320/include" -O2 -Wall -std=gnu99 -mtune=core2 -c C_R_HoltWinters.c -o C_R_HoltWinters.o gcc -m64 -shared -s -static-libgcc -o C_R_HoltWinters.dll tmp.def C_R_HoltWinters.o -Ld:/RCompile/r-compiling/local/local320/lib/x64 -Ld:/RCompile/r-compiling/local/local320/lib -LC:/PROGRA~1/R/R-32~1.2/bin/x64 -lR
Load it in R:
dyn.load('C_R_HoltWinters.dll')
Check availability
> getLoadedDLLs() Filename Dynamic.Lookup base base FALSE utils C:/Program Files/RRO/R-3.2.2/library/utils/libs/x64/utils.dll FALSE methods C:/Program Files/RRO/R-3.2.2/library/methods/libs/x64/methods.dll FALSE RevoUtilsMath C:/Program Files/RRO/R-3.2.2/library/RevoUtilsMath/libs/x64/RevoUtilsMath.dll TRUE grDevices C:/Program Files/RRO/R-3.2.2/library/grDevices/libs/x64/grDevices.dll FALSE graphics C:/Program Files/RRO/R-3.2.2/library/graphics/libs/x64/graphics.dll FALSE stats C:/Program Files/RRO/R-3.2.2/library/stats/libs/x64/stats.dll FALSE tools C:/Program Files/RRO/R-3.2.2/library/tools/libs/x64/tools.dll FALSE internet C:/PROGRA~1/RRO/R-32~1.2/modules/x64/internet.dll TRUE (embedding) (embedding) FALSE C_R_HoltWinters C:/scripts/R/C_R_HoltWinters.dll TRUE
Oh yes, there it is. So, just for poops and chuckles I run a test on it:
> is.loaded('C_R_HoltWinters') [1] FALSE > is.loaded("C_R_HoltWinters") [1] FALSE > is.loaded(C_R_HoltWinters) Error in is.loaded(C_R_HoltWinters) : object 'C_R_HoltWinters' not found
Okay ... he should be there, but it is not. Maybe he knows something that I am not doing, so I am trying to run it:
> .C("C_R_HoltWinters", blahblahblah) Error in .C("C_R_HoltWinters") : C symbol name "C_R_HoltWinters" not in load table > .Call("C_R_HoltWinters", blahblahblah) Error in .Call("C_R_HoltWinters") : C symbol name "C_R_HoltWinters" not in load table
But when I download another c-code called foo and run it, it works fine.
Why can't R refer to C_R_HoltWinters.dll? Will it also break if I put it in a bag?
thanks