ARIMA Model Residuals correlation in small sample? I'm working with simulated data, trying to estimate the best ARIMA model for a time series consisting of 100 observations.
This is the original data. 

First of all, I performed a combination of ADF and KPSS (following what was posted here and here) tests which suggested that the original series was not stationary and had a unit root.

After differentiating once, I ran the tests again and no indication of unit root was shown. In addition to that, I couldn't reject the null for the KPSS test.

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!
 A: Ok I did model your series.


*

*Unit root test, i did a adf test and chose the trend equation since your serie seems to have one, i didnt reject h0 so unit root, I diff the serie and do an other adf test that tells me the serie is now stationnary.




I then used minimization criterion method and it chooses a MA(3) but on diagnostic test for residual, I had autocorrelations so in this case you rise order of p and/or q until you don't have anymore.
It led me to an ARMA(1,5) on the diff serie (so an ARIMA(1,1,5) on original serie)
Here the estimates with all coefficients significant at 10%:

And here the diagnostic test for residuals, LB test until the 24th order to confirm no autocorrelation at 5%:


I have a fairly normal distribution that can allow you to interpret your confidence interval.
For your information, I used SAS 9.4 as you can see and PROC ARIMA.
The methodology is always the same as Box-Jenkins: Statio the serie, find the optimal order based on criterion of your choice (AIC/BIC for explanatory and MSE/others for predictive) and diagnostic test on residuals, normality hypothesis isnt always mandatory if you don't interpret confidence intervals
A: I took your 100 values into AUTOBOX a piece of software that I helped to develop and it suggested a model of (2,0,0) . The residual plot looked clean  with a satisfactory acf of the residuals here  siggesting sufficiency.
The tests that your are using are not always as bad as you present them to be here BUT on the other hand I don't find any of them to be of any substantive value when you have data that has pulses, level shifts, seasonal pulses , local time trends , changes in parameters over tine , changes in error variance over time or any other Gaussian Violation .
Your data is adequately modeled with an ar(2) process WITHOUT any differencing... anything else is unwarranted and superfluous. Can you disclose what the parameters you specified for your simulation ?

So this is the first question: is it possible that this effect is due to the small sample, disregardless of the model chosen?
I have no idea ...
Or is it possible that I need to differentiate again in order to eliminate this issue?
I don't think u need an other order of differencing
