This files allows to estimate and predict using the structural Theta method (see Sbrana and Silvestrini (IJF, 2024)).
The prediction provides both point forecasts and prediction intervals.
The code takes as input the data, the steps and the frequency and returns, as output, a list with the point forecasts (the expected values) and the prediction intervals for the percentages (85%,95%).
If you have any question please email to : giacomo.sbrana@neoma-bs.fr
THETAstructural<-function(y,steps,freq){
if(freq!=1){
s=freq;
for(t in 1:length(y)){if(y[t]!=0){y<-y[t:length(y)];break}}
w<-rep(1/(2*s),s+1);w[2:s]<-1/s
cma<-matrix(NA,length(y),1);
for(g in 1:(length(y)-s)){cma[g+s/2]<-sum(w*y[g:(g+s)])};
residuals<-y/cma
sfactors<-c();for(seas in 1:s){
sfactors[seas]<-mean(na.omit(residuals[seq(seas,length(y)-s+seas,by=s)]))}
sfactors<-sfactors*s/sum(sfactors)
sfactout<-rep(sfactors,length(y)+steps)[(length(y)+1):(length(y)+steps)]
y<-y/rep(sfactors,ceiling(length(y)/s))[1:length(y)]}
likelihood<-function(parameters){q=abs(parameters[1]);co=abs(parameters[2])
state<-p<-rep(0,length(y)+1)
state[1]=y[1]
p[1]=10000
k<-inn<-rep(0,length(y))
k[1]=p[1]/(p[1]+1);
sigmae=0
p[2]=p[1]-p[1]*k[1]+q;
state[2]=co+state[1]+k[1]*inn[1]
for(t in (2):length(y)){
k[t]=p[t]/(p[t]+1);
inn[t]=y[t]-state[t]
sigmae=sigmae+inn[t]^2/(p[t]+1)
p[t+1]=p[t]-p[t]*k[t]+q;
state[t+1]=co+state[t]+k[t]*inn[t]
}
sum(log(p+1))+length(y)*log(sigmae/length(y))}
results<-optim(c(.01,1),likelihood)
q=abs(results[[1]][1]);co=abs(results[[1]][2])
state<-p<-rep(0,length(y)+1)
state[1]=y[1]
p[1]=10000
k<-inn<-rep(0,length(y))
k[1]=p[1]/(p[1]+1);
sigmae=0
p[2]=p[1]-p[1]*k[1]+q;
state[2]=co+state[1]+k[1]*inn[1]
for(t in (2):length(y)){
k[t]=p[t]/(p[t]+1);
inn[t]=y[t]-state[t]
sigmae=sigmae+inn[t]^2/(p[t]+1)
p[t+1]=p[t]-p[t]*k[t]+q;
state[t+1]=co+state[t]+k[t]*inn[t]
}
PointF<-rep(0,steps);PointF[1]=state[length(y)+1];
for( t in 2:steps){PointF[t]=co*(t-1)+state[length(y)+1]}
if(freq!=1){PointF=PointF*sfactout}
H<-sigmae/(length(y)-1);P<-tail(p,1)*H;Eta<-q*H;
Inter<-c();Inter[1]=P;
for(j in 2:steps){Inter[j]=Inter[j-1]+Eta};
Interv<-c();
for(j in 1:steps){Interv[j]=Inter[j]+H};
prob1=.85;prob2=.95;
lower1<-PointF-qnorm((1+prob1)/2)*sqrt(Interv);
for(d in 1:steps){if(lower1[d]<0){lower1[d]=0}};
upper1<-PointF+qnorm((1+prob1)/2)*sqrt(Interv);
lower2<-PointF-qnorm((1+prob2)/2)*sqrt(Interv);
for(d in 1:steps){if(lower2[d]<0){lower2[d]=0}};
upper2<-PointF+qnorm((1+prob2)/2)*sqrt(Interv);
list(mean=PointF,lower=cbind(lower1,lower2),upper=cbind(upper1,upper2))
}