#Recursive formula ( Equation (2.6) in 'Estimated Wold Representation and Spectral Density-Driven Bootstrap for Time Series') 
#Input: Vector of the Fourier coefficients of log f, f a spectral density estimator. In the notation of the paper (a_0,...,a_m)
#Output: A list which contains the estimated Wold coefficients (c_0,...,c_m) and sigma.
#This function computes the resulting Wold coefficients by using the Fourier coefficients of a logaritzmized spectral density estimator. 
#The coefficients are computed up to the order of the input. To compute (c_0,...,c_M) only (a_0,...,a_M) is necessary.
cknfn=function(akn)
{
  tmp=double(length(akn))
  tmp[1]=1
  for(j in 1:(length(akn)-1))
  {
    tmp[j+1]=sum((1-(0:(j-1))/(j))*akn[j-(0:(j-1))+1]*tmp[1:j])
  }
  return(list(ckn=tmp,sigma=exp(akn[1]/2)*sqrt(2*pi)))
}

#Recursive formula (2.7 in 'Estimated Wold Representation and Spectral Density-Driven Bootstrap for Time Series') 
#Input: Vector of the Fourier coefficients of log f, f a spectral density estimator. In the notation of the paper (a_0,...,a_m)
#Output: A list which contains the estimated coefficients (b_0,...,b_m) of the autoregressive representation and sigma.
#This function computes the resulting coefficients of the AR representation by using the Fourier coefficients of a logaritzmized spectral density estimator.
#The coefficients are computed up to the order of the input. To compute (b_0,...,b_M) only (a_0,...,a_M) is necessary.
bknfn=function(akn)
{
  tmp=double(length(akn))
  tmp[1]=-1
  for(j in 1:(length(akn)-1))
  {
    tmp[j+1]=-sum((1-(0:(j-1))/(j))*akn[j-(0:(j-1))+1]*tmp[1:j])
  }
  return(list(bkn=tmp,sigma=exp(akn[1]/2)*sqrt(2*pi)))
}




###Example
require(stats)
require(iosmooth) #Some spectral density estimator
require(datasets)
x=datasets::LakeHuron
n=length(x)
Fouriern=2^(ceiling(log(n,base = 2))+4)
FourierFrequencies=(0:(Fouriern-1))/Fouriern*(2*pi)

fhat=pmax(1/Fouriern,iosmooth::iospecden(x,kernel=c('SupSm'),x.points=FourierFrequencies)$y) #Spectral density estimation
#Be careful, some spectral density estimators return negative values. 

akn=Re(fft(log(fhat))/Fouriern) #Computing Fourier coefficients of log f given log f at Fourier frequencies 0 to 2pi.
ckn=cknfn(akn[1:ceiling(sqrt(length(akn))+1)]) #Computing moving average coefficients

arima.sim(model = list(ma=ckn$ckn[-1]),n = n,rand.gen = rnorm,sd=ckn$sigma) #Generating a time series with spectral density fhat using
#Gaussian noise and the moving average representation


bkn=bknfn(akn[1:ceiling(sqrt(length(akn))+1)]) #Computing autoregressive coefficients

arima.sim(model = list(ar=bkn$bkn[-1]),n = n,rand.gen = rnorm,sd=ckn$sigma) #Generating a time series with spectral density fhat using
#Gaussian noise and the autoregressive representation