Fourier feature for multi seasonality I am trying to forecast a univariate time series with multiple seasonality

*

*Something like this:
library(forecast)
fit <- auto.arima(y, seasonal=FALSE, xreg=fourier(y, K=8))


*Based on visualization my data show frequency of 7, 31 and 365. It is answered in few stack exchange answers that The value of K for fourier term can be chosen by minimizing the AIC.
How would the pseudo code look like to find K for for each frequency ( 7, 31 and 365), to minimize AIC.

You can simply iterate over reasonable numbers of Fourier orders and extract the AIC values, like this:
library(forecast)
taylor_ts <- ts(as.vector(taylor), frequency=336)

max_k <- 10
AICs <- rep(NA,max_k+1)
AICs[1] <- auto.arima(taylor_ts, seasonal=FALSE)$aic

pb <- winProgressBar(max=max_k)
for (kk in 1:max_k) {
    setWinProgressBar(pb,kk,paste(kk,"of",max_k))
    model <- auto.arima(taylor_ts, seasonal=FALSE, xreg=fourier(taylor_ts, K=kk))
    AICs[kk+1] <- model$aic
}
close(pb)

plot(0:max_k,AICs,type="o", pch=19, las=1)


Alternatively, consider defining your time series as an msts object and using a method that is explicitly built for forecasting series with multiple seasonalities.
