Time Series: Confused about identification of (possibly?) an ARMA(p,q) model this is my first ever question on a website i use frequently!
This time series has given me much trouble over the last couple of days even after extensive googling, I suppose with TS theres no two series that look the same, hence my confusion! 
I have the following ACF and PACF plots for my data:

After modelling NUMEROUS models in R (ARMA(p,q), with all combinations of p=0,1,2,3,4 and q=0,1,2,3,4) and looking at their theoretical acf plots, it seems R's auto.arima gave the best model with p=q=2. 
But in my work, i don't want to just say i used R and inferred nothing from the ACF plots haha!
I understand that the slightly significant PACF at lag 1 indicates an AR component(?) and that as both plots 'slope off' it suggests a ARMA(p,q) model, but is there anything else im missing????
Also; the ts (in my eyes) seems seasonal, however ironically, no amount of differencing makes a difference, the plot still looks the same no matter what d equals. (R also tells me its not seasonal)
In short:


*

*What does the ACF plots tell me (main)

*Am i right in not differencing the data


Thanks a bunch!!
 A: The ACF plot shows a seasonality with a period of about 12 (note the sinusoidal pattern). So the first plot to do is a seasonplot.
seasonplot(ts(foo,frequency=12))


However, that does not look very promising.
So, let's look at your time series plot once more. We do see the seasonality, but we also see that the peaks alternate between being high and low. So there does not only seem to be a yearly seasonality, but one of length 2 years as well. So let's look at a seasonplot for a period of 24 months:
seasonplot(ts(foo,frequency=12))


This looks much more regular than the one above. So let's work with a frequency argument of 24. If we now fit a model using auto.arima(), R indeed gives us a seasonal ARIMA model:
model <- auto.arima(ts(foo,frequency=24))
plot(forecast(model,h=24))


I would rather trust auto.arima() to create a forecast based on minimizing AIC, than trying the older Box-Jenkins approach of trying to parse (P)ACF plots. Note that forecast() gives you nice prediction intervals. Take a look at ?forecast.Arima (note the capitalization). Finally, note that in-sample fit is not a good guide to out-of-sample accuracy. Instead, consider using a holdout sample.
A: Trying a list of arima which assume that there are no deterministic structure is frequently inadequate as the correlogram should be calculated from residuals using a model that controls for intervention administration, otherwise the intervention effects are taken to be Gaussian noise, underestimating the actual autoregressive effect. See @AdamO's reflections here @Interrupted Time Series Analysis - ARIMAX for High Frequency Biological Data?.
In a similar ( but different way ) trying to parse/match the original acf/ccf to a set of models (using mid 60's procedures) is also flawed for the same reason ( untreated deterministic structure)
This explains why the following paradigm works well because it integrates both ARIMA (memory) and latent deterministic structure while validating homogeneity/constancy of model parameters and model error variance over time.
https://autobox.com/pdfs/ARIMA%20FLOW%20CHART.pdf
