Modelling moving holiday effects in forecasting I have researched multiple related questions(here, here) but it lacks detailed context and solutions.
My goal is to improve my daily sales forecast accuracy after having incorporated a simple holiday dummy for lunar new year. 
y <- msts(train$Sales, seasonal.periods=c(7,365.25))
# precomputed optimal fourier terms
bestfit$i <- 3
bestfit$j <- 20
z <- fourier(y, K=c(bestfit$i, bestfit$j))
fit <- auto.arima(y, xreg=cbind(z,train_df$cny), seasonal=FALSE)
# forecasting
horizon <- length(test_ts)
zf <- fourier(y, K=c(bestfit$i, bestfit$j), h=horizon)
fc <- forecast(bestfit, xreg=cbind(zf,test_df$cny), h=horizon)
plot(fc, include=365, type="l", xlab="Days", ylab="Sales", main="Comparing arimax forecast and actuals")
lines(test_ts, col='green')


However, this does not reflect the lagged effect of the holiday.

An approach will be to model the effects with a continuous variable(fitted to the effect curve above), but will like to heard other suggestions.
 A: Some points that might help you. 1) Fourier approaches inject structure into the forecasts that might not be present in the data but induced by the equation. You should examine the model residuals carefully 2) One needs to detect the optimal lead and lag around each possible holiday/event and develop a model that includes these effects. As usual there may also be monthly effects/level shift effects/multiple trends as well as anomalous observations. Additionally certain days of the month may be special and even week-in-the-month may play a role. I have implemented a lot of this pre-testing into a commercial package called AUTOBOX which was cited by @forecaster in one of your citataions.
A: You could 


*

*include a number of dummies for the dates around the lunar new year; or

*do as you said (fit a curve).


The first option allows for a lot of flexibility, but that might lead to overfitting. If you have a low number of years in your sample, 1. would probably introduce quite a lot of noise extra to reflecting any signal. It could be nice to apply some form of regularization (shrinkage) to prevent overfitting, but that would be technically challenging.
The second option is less flexible and you need to decide on the functional form, but it is less prone to overfitting.
Thansfer function modelling could probably be a good idea, but I cannot say much about it from my own experience.
