En este programa haremos predicciones con el FOU 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]).
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)

###############################################################################################################
 

kk=1
xfou=rep(0,m)
 

for (ii in 1:m)
{


xx=x[1:(length(x)-1-m+ii)]
  
  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)
 
 lg=(2*mean(x^2)/(og^2*gamma(2*Hg+1)))^(-1/(2*Hg))#lambda gorro

n=length(xx)


  
  
  
     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)
  
      }
     
      og^2*Hg*f(lg*t)/(2*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
  
  
  xfou[ii]=pred
  print(ii) 
}


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

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

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