En este programa haremos predicciones con el FOU(l1,l1) para una serie de datos estacionaria y centrada.
Las predicciones las haremos a un paso. Estimaremos el modelo hasta cada dato y luego predecimos 
el siguiente. 


Ingresar x (la serie debe ser centrada y estacionaria).
Ingresar m (cantidad de valores a predecir).
Ingresar Te (intervalo [0,T]).
Ingresar aa y bb (parámetros de la función w, aa es el a del teorema, aa+bb es el b del teorema).
Las estimaciones de los parámetros se obtienen corriendo el programa para cualquier valor de m
y se los llama mediante Hg, og y lg (H gorro, sigma gorro y lambdagorro)

###############################################################################################################
 
aa=4
bb=3
kk=1
xfou2_l1_l1=rep(0,m)#acá pondremos lakegorro (las predicciones a un paso)
 

for (ii in 1:m)
{


xx=x[1:(length(x)-1-m+ii)]
  #Ahora ajusto los parámetros del FOU(2)
  
  n=length(xx)
  
  a=c(0.4829629131445341,-0.8365163037378077,0.2241438680420134,0.1294095225512603)/sqrt(2)
  #a=c(-1,2,-1)/4#dio Hg negativo
  #a=c(-1,3,-3,1)/8#dio Hg negativo
  #a=c(1,-4,6,-4,1)/16#dio Hg negativo
  #a=c(-1,5,-10,10,-5,1)/32#dio Hg negativo
  #a=c(-1,6,-15,20,-15,6,-1)/64#dio Hg negativo
  #a=c(-1,7,-21,35,-35,21,-7,1)/128#dio Hg negativo
#  a=c(-1,8,-28,56,-70,56,-28,8,-1)/256#dio Hg negativo

  
 
  p=length(a)
  del=Te/n
  
  v1=rep(0,(length(xx)-p+1))
  v2=rep(0,(length(xx)-2*p))
  #Ahora calculamos el filtro de longitud 2p
  a2=rep(0,p)
  a2=rbind(a,a2)
  a2=a2[1:(2*p-1)]
  
  
  for(i in 1:(length(xx)-p+1))
  {
    v1[i]=(sum(a*xx[i:(i+p-1)]))^2
  }
  va=mean(v1)
  
  for(i in 1:(length(xx)-2*p+2))
  {
    v2[i]=(sum(a2*xx[i:(i+2*p-2)]))^2
  }
  va2=mean(v2)
  
  Hg=log2(va2/va)/2# Hg=Hgorro
  b=rep(0,length(a))
  
  for(j in 1:length(a))
  {
    A=seq(1,p)
    A=abs(A-j)^(2*Hg)
    b[j]=sum(a*A)
  }
  
  og=(-2*va/(del^(2*Hg)*sum(a*b)))^0.5# sigma gorro (ogorro)


n=length(xx)

  U_Tn=function(l)
  {
    j=seq(Te/n,Te,Te/n)
    w=abs(j)^aa/(1+abs(j)^(aa+bb))
    esp_dens=og^2*gamma(2*Hg+1)*sin(Hg*pi)*abs(j)^(3-2*Hg)/(2*pi*(l^2+j^2)^2)    
    
    Idelta=rep(NA,n)
    for (i in 1:n)
    {
      Idelta[i]=Te/(2*pi)*(mean(xx*cos(j*i*Te/n))^2+mean(xx*sin(j*i*Te/n))^2)
    }
    (Te/(2*pi*n))*sum((log(esp_dens)+Idelta/esp_dens)*w)
    
  }
  
  #lg=optim(c(0.1),U_Tn)$par
  lg=optim(c(0.1),U_Tn, lower=0.01,upper=1.5,method ="L-BFGS-B")$par
 
  
  
     Covar=function(t)
    {
      f=function(x)
      {
      A=rep(0,170)
      B=rep(0,170)
      for(i in 1:170)
       {
        A[i]=x^(i-1)/(factorial(i-1)*(i-1+2*Hg))
        B[i]=(-x)^(i-1)/(factorial(i-1)*(i-1+2*Hg))
       }
       gamma(2*Hg)*(exp(x)+exp(-x))-exp(x)*x^(2*Hg)*sum(B)-exp(-x)*x^(2*Hg)*sum(A)
  
      }
      fp=function(x)
      {
            A=rep(0,170)
      for(i in 1:170)
       {
        A[i]=x^(i-1)/(factorial(i-1)*(i-1+2*Hg))
       }
      f(x)-2*exp(-x)*gamma(2*Hg)+2*exp(-x)*x^(2*Hg)*sum(A)-2*x^(2*Hg-1)
  
      }
      
      og^2*Hg*((2-2*Hg)*f(lg*t)+lg*t*fp(lg*t))/(4*lg^(2*Hg))
      
    }
    
  
  
  n=length(xx)
  
  covarvector=c(rep(0,(n+kk)))
  for(i in 2:(n+kk))
  {
    
    
    covarvector[i]=Covar((i-1)*Te/(n+kk-1))
  }
  covarvector[1]=og^2*Hg*(1-Hg)*gamma(2*Hg)/(2*lg^(2*Hg))#acá ponemos la mitad de la varianza ya que al sumar t(C) debe quedar en la diagonal la varianza)
  
  C=matrix(rep(0),(n+kk),(n+kk))
  C[,1]=covarvector[1:(n+kk)]
  
  
  
  for(j in 1:(n+kk-1))
  {
    C[,(j+1)]=c(rep(0,j),covarvector[1:(n+kk-j)])
   # print(j)
  }
  
  
  C=C+t(C)
  
  #sum(abs(C-t(C)))#con esto controlo que la matriz haya quedado simétrica
  
  #min(eigen(C)$values)#así verifico que es definida positiva
  #sum(eigen(C)$values<0)
  pred=C[(n+1):(n+kk),1:n]%*%solve(C[1:n,1:n])%*%xx
  
  
  xfou2_l1_l1[ii]=pred
  print(ii) 
}


plot(x[(length(x)+1-m):length(x)],type="o")
lines(xfou2_l1_l1,col="red",type="o")

#Agregamos ahora el indice de Willmott de las predicciones
obs=x[(length(x)+1-m):length(x)]
pre=xfou2_l1_l1

W2=1-sum((pre-obs)^2)/sum((abs(pre-mean(obs))+abs(obs-mean(obs)))^2)
RMSE=sqrt(mean((pre-obs)^2))

W1=1-sum(abs((pre-obs)))/sum((abs(pre-mean(obs))+abs(obs-mean(obs))))
MAE=mean(abs(pre-obs))

W2
RMSE
W1
MAE