Результат оптимизации случайный

Question

У меня есть скрипт, который запускает ARIMA, устанавливая весы на ошибки. Скрипт работает отлично, но каждый раз он запускается, даже используя одну и ту же серию, он выводит разные прогнозы. Я просмотрел весь код и не могу найти, где проблема. Я был бы очень признателен, если бы кто-то мог быстро взглянуть и указать, где я ошибся.

M<-matrix(c("08Q1", "08Q2", "08Q3", "08Q4", "09Q1", "09Q2", "09Q3", "09Q4", "10Q1", "10Q2", "10Q3", "10Q4", "11Q1", "11Q2", "11Q3", "11Q4", "12Q1", "12Q2", "12Q3", "12Q4", "13Q1", "13Q2", "13Q3", "13Q4", "14Q1", "14Q2", 79160.56, 91759.73, 91186.48, 106353.82, 70346.47, 80279.15, 82611.60, 131392.72, 93798.99, 105944.78, 103913.13, 154530.69, 110157.40, 117416.09, 127423.42, 156752.00,120097.81, 121307.75, 115021.12, 150657.83, 113711.53, 115353.14, 112701.98, 154319.18,116803.54, 118352.54),ncol=2,byrow=FALSE) 
deltaT<-1/4 
horiz<-4 
startY<-c(8,1) 
aslog<-"y" 
Nu<-M[,length(M[1,])] 
Nu<-as.numeric(Nu) 
Nu<-ts(Nu,deltat=deltaT,start=startY) 
Mdates<-as.character(M[,1]) 
if(aslog=="y") 
    {N<-log(Nu)} else 
    {N<-Nu} 
library(forecast) 
library(tseries) 
max.sdiff <- 3 
arima.force.seasonality <- "n" 
    fweight <- function(x) 
    { 
    PatX <- 0.5+x 
    return(PatX) 
    } 
    integ1 <- integrate(fweight, lower = 0.00, upper = 1) 
    valinteg <- 2*integ1$value 
integvals <- rep(0, length.out = length(N)) 
for (i in 1:length(N)) 
{ 
integi <- integrate(fweight, lower = (i-1)/length(N), upper= i/length(N)) 
integvals[i] <- 2*integi$value 
} 

kpssW <- kpss.test(N, null="Level") 
ppW <- tryCatch({ppW <- pp.test(N, alternative = "stationary")}, error = function(ppW) {ppW <- list(error = "TRUE", p.value = 0.99)}) 
adfW <- adf.test(N, alternative = "stationary", k = trunc((length(N)-1)^(1/3))) 
if(kpssW$p.value < 0.05 | ppW$p.value > 0.05 | adfW$p.value > 0.05) {ndiffsW = 1} else {ndiffsW = 0} 
aaw <- auto.arima(N, max.D= max.sdiff, d=ndiffsW, seasonal=TRUE, 
    allowdrift=FALSE, stepwise=FALSE, trace=TRUE, seasonal.test="ch") 
orderWA <- c(aaw$arma[1], aaw$arma[6] , aaw$arma[2]) 
orderWS <- c(aaw$arma[3], aaw$arma[7] , aaw$arma[4]) 
if(sum(aaw$arma[1:2])==0) {orderWA[1] <- 1} else {NULL} 
if(arima.force.seasonality == "y") {if(sum(aaw$arma[3:4])==0) {orderWS[1] <- 1} else {NULL}} else {NULL} 
stAW <- Arima(N, order= orderWA, seasonal=list(order=orderWS), method="ML") 
parSW <- stAW$coef 
WMAEOPT <- function(parSW) 
{ 
    ArimaW <- Arima(N, order = orderWA, seasonal=list(order=orderWS), 
    include.drift=FALSE, method = "ML", fixed = c(parSW)) 
    errAR <- c(abs(resid(ArimaW))) 
    WMAE <- t(errAR) %*% integvals 
    return(WMAE) 
} 
OPTWMAE <- optim(parSW, WMAEOPT, method="SANN", control = list(fnscale= 1, maxit = 5000)) 
parS3 <- OPTWMAE$par 
ArimaW1 <- Arima(N, order = orderWA, seasonal=list(order=orderWS), 
    include.drift=FALSE, method = "ML", fixed = c(parS3)) 
fArimaW1 <- forecast(ArimaW1, h=8, simulate= TRUE, fan=TRUE) 
if (aslog == "y") {fArimaWF <- exp(fArimaW1$mean[1:horiz])} else {fArimaWF <- fArimaW1$mean[1:horiz]} 
plot(fArimaW1, main = "ARIMA Forecast", sub="blue=fitted, red=actual") # ylim=c(17, 20) 
    lines(N, col="red", lwd=2) 
    lines(ts(append(fitted(ArimaW1), fArimaW1$mean[1]), deltat=deltaT, start = startY), 
    col= "blue", lwd = 2) # makes the graph look nicer 
if (aslog == "y") {ArimaALT <- exp(fArimaW1$mean[1:horiz])} else {ArimaALT <- fArimaW1$mean[1:horiz]} 
start(fArimaW1$mean) -> startF 
ArimaALTf <- ts(prettyNum(ArimaALT, big.interval = 3L, big.mark = ","), deltat = deltaT , start= startF) 
View(ArimaALTf, title = "ARIMA forecast") 
summary(ArimaW1) 

Редактировать

я только что нашел, где это идет не так. Но я не понимаю, почему.

OPTWMAE <- optim(parSW, WMAEOPT, method="SANN", control = list(fnscale= 1, maxit = 5000)) 

Это где он дает различные значения

Благодарим за Ваше время

Answer 1

От help("optim") (курсив мной):

Метод "Санна" по умолчанию а вариант моделируемого отжига, приведенный в Belisle (1992). Имитированный отжиг относится к классу стохастик глобальные методы оптимизации.

Используйте set.seed, чтобы получить воспроизводимые результаты.