After that, I used auto.arima (focusing on BIC criteria for model selection). I also created a pool of models and ranked them according to BIC (which, as expected, did not match the auto.arima suggestion). This is a data.frame that I created with each model's specification and it's corresponding BIC value.
After that, I started testing the residuals of each model: Ran several normality tests (JB, SW, K-S) as well as plots for visual identification of 'normality'. So far, no issues here for any model.
The problem appears when I run the Ljung-Box test to test for correlation of the residuals.
For every model that I tested, the residuals appear to be correlated. For one particular model, I even tested for the presence of ARCH effect (using Q and LM tests) and there were signs of such effect.
So this is the first question: is it possible that this effect is due to the small sample, disregardless of the model chosen? Or is it possible that I need to differentiate again in order to eliminate this issue? Among the models in the table, there are some that require differentiation a second time. Did that, tested for autocorrelation (again, using Ljung-Box test) in the residuals and couldn't reject the null. All this is shown in the next image:
Am I failing to identify the presence of unit root after differentiation for the first time? Is it auto.arima as well as my function that created a pool of models failing to identify the need for a second differentiation? Or is there any explanation given the size of the data for this issue with the autocorrelation of residuals?
If necessary, data can be found here.
Thanks in advance!