I have weekly data of yearly seasonality around 3.5 years. By using default seasonal test ocsb or ch test in python pmdarima it is not able to give right D differencing which is resulting in higher errors. I am using exogenous 0/1 dummy categories of major peaks of new year sales ,etc etc .
Can I use the R forecast library seasonal strength formula to determine the D Order .I am doing regression with best suitable dummy 0/1 categories and passing the residuals to STL decomposition and applying R forecast default seasonal test strength formula score (.64 cutoff) to find the right D (using 0 or 1 only) and using that D in auto.arima along with those best suitable dummy categories exogenous. Does this make sense ?