> kt
 [1] 16.63292944 17.10746185 16.03326222 16.24631736 14.67784276
 [6] 14.54126450 14.12032640 14.08797796 13.69164459 13.32990466
[11] 13.15672641 12.57994701 12.50856702 12.49292845 11.91980041
[16] 11.84854379 11.69906707 11.27068958 11.46658995 11.07480840
[21] 11.71810389 10.24234202  9.98624154  9.69563628  8.77524312
[26]  8.42607407  8.03839654  7.46459155  6.90364116  6.46802360
[31]  6.18895143  5.66047744  5.06493622  4.87223437  4.47053611
[36]  4.31558216  3.80968168  3.54183060  3.32698878  3.01142929
[41]  2.61087797  2.16596582  1.84685422  1.73129589  1.44568725
[46]  1.38868046  1.05136432  0.65361666  0.31751012 -0.06360237
[51] -0.56454474 -0.78109888 -0.82625274 -0.81258738 -1.01088682
[56] -0.83005175 -1.05284966 -0.88126429 -0.88002998 -0.73160603
[61] -0.73800657 -0.82408653 -0.84909859 -0.86868214 -0.82202275
[66] -0.67518149

> ktF.fit <- auto.arima(kt, d=1, trace=T)

ktF.for <- forecast(ktF.fit,h=horizonte,level=c(95))

set.seed(1)
ktF.coef <- rmvnorm(nsimul, ktF.fit$coef, ktF.fit$var.coef, method="chol")

set.seed(1)
NkF <- (horizonte+2)*nsimul
normalkF <- matrix(rnorm(NkF,0,sqrt(ktF.fit$sigma2)),(horizonte+2),nsimul)

ktF <- matrix(NA,(horizonte+4),nsimul)

ktF[1,] <- ktF.fit$x[tiempo-2]
ktF[2,] <- ktF.fit$x[tiempo-1]
ktF[3,] <- ktF.fit$x[tiempo]

for(i in 4:(horizonte+4)){
  for(j in 1:nsimul){
    ktF[i,j] <- ktF[i-1,j] +
      ktF.coef[j,4]*(1-ktF.coef[j,1]-ktF.coef[j,2]) + #drift
      ktF.coef[j,1]*(ktF[i-1,j]-ktF[i-2,j]) + #Ar(1)
      ktF.coef[j,2]*(ktF[i-2,j]-ktF[i-3,j]) + #Ar(2)
      normalkF[i-2,j] +
      ktF.coef[j,3]*normalkF[i-3,j] #MA(1)
  }
}

ktF1<-ktF[4:(horizonte+4),1:nsimul]

ic.ktF<-matrix(NA,(horizonte+1),3)
for(i in 1:(horizonte+1)){
  ic.ktF[i,]<-quantile(ktF1[i,],probs=c(0.025,0.5,0.975),
                       type = 8,names=FALSE,na.rm=T)
}