lower <- c()
upper <- c()
c_j <- alpha+beta*seq(1,5)
v <- 0

for (i in seq(1, 6)){
  if(i == 1){lower = c(lower,
              qnorm(0.025, mean = l_t[length(l_t)] + i*b_t[length(b_t)],
              sd = sqrt(sigma2)))
  upper = c(upper,
            qnorm(0.975, mean = l_t[length(l_t)] + i*b_t[length(b_t)],
                  sd = sqrt(sigma2)))
  }
  else{

  v = sigma2*(1 + sum(c_j[1:i-1]^2))

  lower = c(lower,
            qnorm(0.025, mean = l_t[length(l_t)] + i*b_t[length(b_t)],
                  sd = sqrt(v)))
  upper = c(upper,
            qnorm(0.975, mean = l_t[length(l_t)] + i*b_t[length(b_t)]
                  sd = sqrt(v)))
  }
}

forecast_analytic1 <- cbind(prediction, lower, upper)
colnames(forecast_analytic1) <- c('Forecast', 'Lower', 'Upper')`