ARMA modeling in R I have a time series. I want to model it using ARMA, which will be used for forcasting. 
In R I am using arima() function to get the coefficients. But arima() requires order(p,d,q) as input. What is the simplest way in R to arrive at a good value for p and q (with d = 0) so that I don't overfit?
 A: Simplest way to arrive at values for $p$ and $q$ is using auto.arima function from package forecast. There is no simplest way in any statistical package to arrive at good values. The main reason for that is that there is no universal definition of good. 
Since you mention overfitting, one possible way is to fit arima models for different values of $p$ and $q$ and then pick the one which is the best according to your overfitting criteria (out of sample forecasting performance for example). auto.arima does basically the same, you can choose between AIC, AICC and BIC to let auto.arima pick the best model. 
A: One option is to fit a series of ARMA models with combinations of $p$ and $q$ and work with the model that has the best "fit". Here I evaluate "fit" using BIC to attempt to penalise overly complex fits. An example is shown below for the in-built Mauna Loa $\mathrm{CO}_2$ concentration data set
## load the data
data(co2)
## take only data up to end of 1990 - predict for remaining data later
CO2 <- window(co2, end = c(1990, 12))

## Set up the parameter sets over which we want to operate
CO2.pars <- expand.grid(ar = 0:2, diff = 1, ma = 0:2, sar = 0:1,
                        sdiff = 1, sma = 0:1)
## As you are only wanting ARMA, then you would need something like
## pars <- expand.grid(ar = 0:4, diff = 0, ma = 0:4)
## and where you choose the upper and lower limits - here 0 and 4

## A vector to hold the BIC values for each combination of model
CO2.bic <- rep(0, nrow(CO2.pars))

## loop over the combinations, fitting an ARIMA model and recording the BIC
## for that model. Note we use AIC() with extra penalty given by `k`
for (i in seq(along = CO2.bic)) {
    CO2.bic[i] <- AIC(arima(CO2, unlist(CO2.pars[i, 1:3]), 
                            unlist(CO2.pars[i, 4:6])),
                      k = log(length(CO2)))
}

## identify the model with lowest BIC
CO2.pars[which.min(CO2.bic), ]

## Refit the model with lowest BIC
CO2.mod <- arima(CO2, order = c(0, 1, 1), seasonal = c(0, 1, 1))
CO2.mod
## Diagnostics plots
tsdiag(CO2.mod, gof.lag = 36)

## predict for the most recent data
pred <- predict(CO2.mod, n.ahead = 7 * 12)
upr <- pred$pred + (2 * pred$se) ## upper and lower confidence intervals
lwr <- pred$pred - (2 * pred$se) ## approximate 95% pointwise

## plot what we have done
ylim <- range(co2, upr, lwr)
plot(co2, ylab = ylab, main = expression(bold(Mauna ~ Loa ~ CO[2])),
     xlab = "Year", ylim = ylim)
lines(pred$pred, col = "red")
lines(upr, col = "red", lty = 2)
lines(lwr, col = "red", lty = 2)
legend("topleft", legend = c("Observed", "Predicted", "95% CI"),
       col = c("black", "red", "red"), lty = c(1, 1, 2), bty = "n")

