I am trying to build an ARMA GARCH model and since I haven't found much information about the actual process I'll try it here.

So, what have I allready done?

I have tested my data for stationarity and consequently I have differenced it to make it stationary. I have then evaluated several ARIMA specifications and according to AIC,BIC and Q-test I have decided to use AR(1/3,5,6) model, which exhibits almost none autocorrelation among the residuals.

Then, I have checked for ARCH effects using the archlm test and indeed, there were some. Hence, I have examined autocorrelation among the squared residuals from my ARIMA model and there were some significant values (29% at first lag, 13% at the fourth).

Having confirmed the ARCH effect I have decided to include GARCH(1,1) to my previous model and there are two scenarios I have tried.

1. I have constrained the coefficients of the AR terms from the previous estimation to keep the mean-modelling intact as it seeme very promising among the others ARIMA specifications I tried. Although the residuals didn|t have anz signs of autocorrelation, the autocorrealtion of the squared residuals at the first lag decreased only about 0,5% which isn't much of a improvement, is it? Including further lags to the GARCH doesn't help either.
2. I have estimated then new model unconstrained and again only very small improvement and moreover the autocorrelation of the residuals (not squared), which was previously diminished by the AR terms, appears again.

To conclude I am able to identify and fit the ARMA, but I am not able to fit the GARCH. Am I doing it wrong?

PS. One thing which confuses me a bit and which I was not able to find is whether to determine the ARMA model and then use the GARCH model wit appropriate ARMA terms or if I should estimate the ARMA using my dependent variable, then take the residuals of the estimation and use them as dependent variable for the garch modelling?

As @Dan requested here are the plots

My time series

ACF of the time series

ACF of the residuals after the ARMA estimation

ACF of the squared residuals