This files allows to estimate and predict using the double damped trend model (see Sbrana and Yu (JORS, 2024)).
The code contains the estimation of both MSOE and SSOE versions. The prediction part provides both point forecast and prediction intervals.
Both codes take as input the data and the steps and return, 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 : \(\textbf{giacomo.sbrana@neoma-bs.fr}\)
DDTmsoe<-function(y,h){
prob1=.85;prob2=.95;
damped<-matrix(0,length(y),2)
damped[1,1]=y[1]
damped[1,2]=0#y[2]-y[1]
inn<-matrix(0,length(y),1)
fmsoe <- function(param){a<-exp(param[1])/(1+exp(param[1]));b<-exp(param[2])/(1+exp(param[2]));w<-exp(param[3])/(1+exp(param[3]));x<-exp(param[4])/(1+exp(param[4]));
p1=-(1/(4* a^2* b^2))*(1 + b^2 - a^2* (-1 + 3* b^2)* (1 + w) + b^2* x + sqrt(1 + a^4 *(-1 + b^2)^2 *(1 + w)^2 + 2* b^2* (-1 + x) + 8* a* b^3 *x + b^4 *(1 + x)^2 + 2* a^2 *((-1 + b^2)^2* (-1 + w) + b^2* (1 + b^2)* (1 + w)* x)) - sqrt(2)* a^2 *b^2* sqrt((1/(a^4* b^4)*((1 - a^2 *(-1 + 3* b^2) *(1 + w) + b^2* (1 + x))^2 + 2* a^2* b^2* (-1 + w - a^2 *(-1 + b^2)* (1 + w)^2 + b^2* (1 + w)* (1 + x)) + (((-1 + b)^2* (1 + a* (-2 + a + a* w)) + b^2* x)* ((1 + b)^2* (1 + a *(2 + a + a *w)) + b^2 *x) *(1 + a^2* (1 + b^2) *(1 + w) + b^2 *(1 + x)))/sqrt( 1 + a^4 *(-1 + b^2)^2* (1 + w)^2 + 2* b^2* (-1 + x) + 8* a* b^3 *x + b^4 *(1 + x)^2 + 2* a^2 *((-1 + b^2)^2 *(-1 + w) + b^2 *(1 + b^2)* (1 + w) *x)) + 2* (-b^2 + a* (a* (-1 + b^2 *(1 + w - a^2 *(-2 + 3 *b^2)* (1 + w)^2 + b^2* (2 + 3 *w))) + b^3* (2 + 3* a* b *(1 + w)) *x))))))
k1<- p1-w;
k2<- (w+p1*(-1+a^2*(1-p1+w)))/(-1+a*b*(-1+p1-w));
for (t in 2:(length(y))) {
inn[t]<-y[t]-a*damped[t-1,1]-b*damped[t-1,2]
damped[t,1] <- a*damped[t-1,1]+b*damped[t-1,2]+k1*inn[t]
damped[t,2] <- b*damped[t-1,2]+k2*inn[t]
}
if(max(y)>(mean(y[-which.max(y)])+3*sd(y[-which.max(y)]))){inn=inn[-which.max(y)]}
sum(inn^2)}
result<-optim(c(1.4,1,-.5,-1.5),fmsoe)
aa=exp(result[[1]][1])/(1+exp(result[[1]][1]))
bb=exp(result[[1]][2])/(1+exp(result[[1]][2]))
ww=exp(result[[1]][3])/(1+exp(result[[1]][3]))
xx=exp(result[[1]][4])/(1+exp(result[[1]][4]))
p1=-(1/(4* aa^2* bb^2))*(1 + bb^2 - aa^2* (-1 + 3* bb^2)* (1 + ww) + bb^2* xx + sqrt(1 + aa^4 *(-1 + bb^2)^2 *(1 + ww)^2 + 2* bb^2* (-1 + xx) + 8* aa* bb^3 *xx + bb^4 *(1 + xx)^2 + 2* aa^2 *((-1 + bb^2)^2* (-1 + ww) + bb^2* (1 + bb^2)* (1 + ww)* xx)) - sqrt(2)* aa^2 *bb^2* sqrt((1/(aa^4* bb^4)*((1 - aa^2 *(-1 + 3* bb^2) *(1 + ww) + bb^2* (1 + xx))^2 + 2* aa^2* bb^2* (-1 + ww - aa^2 *(-1 + bb^2)* (1 + ww)^2 + bb^2* (1 + ww)* (1 + xx)) + (((-1 + bb)^2* (1 + aa* (-2 + aa + aa* ww)) + bb^2* xx)* ((1 + bb)^2* (1 + aa *(2 + aa + aa *ww)) + bb^2 *xx) *(1 + aa^2* (1 + bb^2) *(1 + ww) + bb^2 *(1 + xx)))/sqrt( 1 + aa^4 *(-1 + bb^2)^2* (1 + ww)^2 + 2* bb^2* (-1 + xx) + 8* aa* bb^3 *xx + bb^4 *(1 + xx)^2 + 2* aa^2 *((-1 + bb^2)^2 *(-1 + ww) + bb^2 *(1 + bb^2)* (1 + ww) *xx)) + 2* (-bb^2 + aa* (aa* (-1 + bb^2 *(1 + ww - aa^2 *(-2 + 3 *bb^2)* (1 + ww)^2 + bb^2* (2 + 3 *ww))) + bb^3* (2 + 3* aa* bb *(1 + ww)) *xx))))))
p2=(ww+p1*(-1 + aa^2 +aa^2* (-p1 + ww)))/(-1-aa*bb+aa*bb*(p1 - ww))
p3=(1+aa^2*p1+(1/(1-p1+ww))-((2*(1+aa*bb+aa^2*p1))/(1+aa*bb+aa*bb*(ww-p1))))/(bb^2)
P=matrix(c(p1,p2,p2,p3),2,2)
F= 1+ aa^2*p1+2*aa*bb*p2+bb^2*p3
k1<- (p1-ww);
k2<- (ww+p1*(-1+aa^2*(1-p1+ww)))/(-1+aa*bb*(-1+p1-ww));
for (t in 2:(length(y))) {
inn[t]<-y[t]-aa*damped[t-1,1]-bb*damped[t-1,2]
damped[t,1] <- aa*damped[t-1,1]+bb*damped[t-1,2]+k1*inn[t]
damped[t,2] <- bb*damped[t-1,2]+k2*inn[t]
}
Tmat<-matrix(bb,2,2);Tmat[1,1]<-aa;Tmat[2,1]<-0;
T=diag(2)
Forec<-c();
for(i in 1:steps){
Forec[i]<-c(aa,bb)%*%T%*%(damped[length(y),]);
T=T%*%Tmat
}
H<-sum(inn^2)/(length(inn)*F);#Eta<-diag(c(ww,xx))*H
T=Tmat;Z=c(aa,bb);Eta=P*H;P=P*H;
Inter<-list();Inter[[1]]=P;for(j in 2:steps){Inter[[j]]=T%*%Inter[[j-1]]%*%t(T)+Eta}
Interv<-c();for(j in 1:steps){Interv[j]=t(Z)%*%Inter[[j]]%*%Z+H};
lower1<-Forec-qnorm((1+prob1)/2)*sqrt(Interv);for(d in 1:steps){if(lower1[d]<0){lower1[d]=0}}
upper1<-Forec+qnorm((1+prob1)/2)*sqrt(Interv);
lower2<-Forec-qnorm((1+prob2)/2)*sqrt(Interv);for(d in 1:steps){if(lower2[d]<0){lower2[d]=0}}
upper2<-Forec+qnorm((1+prob2)/2)*sqrt(Interv)
list(mean=Forec,lower=cbind(lower1,lower2),upper=cbind(upper1,upper2))
}
DDTssoe<-function(y,h){
prob1=.85;prob2=.95;
damped<-matrix(0,length(y),2)
damped[1,1]=y[1]
damped[1,2]=0#y[2]-y[1]
inn<-matrix(0,length(y),1)
fssoe <- function(param){a<-exp(param[1])/(1+exp(param[1]));b<-exp(param[2])/(1+exp(param[2]));k1<-exp(param[3])/(1+exp(param[3]));k2<-exp(param[4])/(1+exp(param[4]));
for (t in 2:(length(y))) {
inn[t]<-y[t]-a*damped[t-1,1]-b*damped[t-1,2]
damped[t,1] <- a*damped[t-1,1]+b*damped[t-1,2]+k1*inn[t]
damped[t,2] <- b*damped[t-1,2]+k2*inn[t]
}
if(max(y)>(mean(y[-which.max(y)])+3*sd(y[-which.max(y)]))){inn=inn[-which.max(y)]}
sum(inn^2)}
result<-optim(c(1.4,1,-.5,-1.5),fssoe)
aa=exp(result[[1]][1])/(1+exp(result[[1]][1]))
bb=exp(result[[1]][2])/(1+exp(result[[1]][2]))
k1=exp(result[[1]][3])/(1+exp(result[[1]][3]))
k2=exp(result[[1]][4])/(1+exp(result[[1]][4]))
for (t in 2:(length(y))) {
inn[t]<-y[t]-aa*damped[t-1,1]-bb*damped[t-1,2]
damped[t,1] <- aa*damped[t-1,1]+bb*damped[t-1,2]+k1*inn[t]
damped[t,2] <- bb*damped[t-1,2]+k2*inn[t]
}
Tmat<-matrix(bb,2,2);Tmat[1,1]<-aa;Tmat[2,1]<-0;
T=diag(2)
Forec<-c();
for(i in 1:steps){
Forec[i]<-c(aa,bb)%*%T%*%(damped[length(y),]);
T=T%*%Tmat
}
H<-sum(inn^2)/length(inn);Eta<-P<-c(k1,k2)%*%t(c(k1,k2))*H
T=Tmat;Z=c(aa,bb);
Inter<-list();Inter[[1]]=P;for(j in 2:steps){Inter[[j]]=T%*%Inter[[j-1]]%*%t(T)+Eta}
Interv<-c();for(j in 1:steps){Interv[j]=t(Z)%*%Inter[[j]]%*%Z+H};
lower1<-Forec-qnorm((1+prob1)/2)*sqrt(Interv);for(d in 1:steps){if(lower1[d]<0){lower1[d]=0}}
upper1<-Forec+qnorm((1+prob1)/2)*sqrt(Interv);
lower2<-Forec-qnorm((1+prob2)/2)*sqrt(Interv);for(d in 1:steps){if(lower2[d]<0){lower2[d]=0}}
upper2<-Forec+qnorm((1+prob2)/2)*sqrt(Interv)
list(mean=Forec,lower=cbind(lower1,lower2),upper=cbind(upper1,upper2))
}