Correlogram in R like in Stata? In STATA I can create a "Correlogram" to find the appropriate lag order in case of time series. E.g.

I know I can use the acf or Acf of the forecast package to calculate the ACF and PACF and to plot it. But how can I get the significance values? I mean the values in the column Prob>Q? Is this implemented in any package/command?
 A: You can calculate the threshold significance for a given coefficient.  To check this, I just did a quick google search.  Found this https://stat.ethz.ch/pipermail/r-help/2009-August/207266.html The code to make the calculation is below.  I hope this helps.
dev.new(height=6, width=4.5)
par(mfrow=c(2,1), mar=c(3,4,0.5,0.5))

testAR1 <- arima.sim(n=300, list(ar=c(0.5)))
test_acf <- acf(testAR1)
test_acf_CritVal <- qnorm((1 + 0.95)/2)/sqrt(test_acf$n.used) #change the 0.95 to the desired CI
abline(h=test_acf_CritVal, col="red") #add the calculated CI to the plot, just to check that it goes in the right spot

test_pacf <- pacf(testAR1)
test_pacf_CritVal <- qnorm((1 + 0.95)/2)/sqrt(test_pacf$n.used) #change the 0.95 to the desired CI
abline(h=test_pacf_CritVal, col="red") #add the calculated CI to the plot, just to check that it goes in the right spot

Now that you have the critical value, you can now compare that to the estimated coefficients for the different lags.
Test_acf_lags <- as.vector(test_acf$lag)
Test_acf_coefs <- as.vector(test_acf$acf)
Test_acf_Signif <- as.integer(Test_acf_coefs>test_acf_CritVal)

Test_pacf_lags <- as.vector(test_pacf$lag)
Test_pacf_coefs <- as.vector(test_pacf$acf)
Test_pacf_Signif <- as.integer(Test_pacf_coefs>test_pacf_CritVal)

R_Correlo_Table_acf <- data.frame("Lag"=Test_acf_lags, "acf"=Test_acf_coefs, "Signif_acf"=Test_acf_Signif)
R_Correlo_Table_pacf <- data.frame("Lag"=Test_pacf_lags, "pacf"=Test_pacf_coefs, "Signif_pacf"=Test_pacf_Signif)
R_Correlo_Table <- merge(R_Correlo_Table_acf, R_Correlo_Table_pacf, all=TRUE)

print(R_Correlo_Table)

EDIT: Added the calculation and plotting of confidence bands for model ID
dev.new(height=6, width=4.5)
par(mfrow=c(2,1), mar=c(3,4,0.5,0.5))

Nsim <- 100
test2AR1 <- arima.sim(n=Nsim, list(ar=c(0.85)))
test2_acf <- acf(test2AR1, ylim=c(-1,1))
Nlag <- length(test2_acf$lag)
Bart_se_acf_frist <- 1/sqrt(Nsim)
Bart_se_acf_rest <- sqrt((1+2*cumsum(test2_acf$acf^2)[-Nlag])/Nsim)
Bart_se_acf <- c(Bart_se_acf_frist, Bart_se_acf_rest)
ConfBands <- qnorm((1-0.05/2)) * Bart_se_acf
lines(test2_acf$lag, ConfBands, col="red") #add the calculated CI to the plot, just to check that it goes in the right spot
lines(test2_acf$lag, ConfBands*-1, col="red")


test2_pacf <- pacf(test2AR1, ylim=c(-1,1))
Nlag <- length(test2_pacf$lag)
Bart_se_pacf_frist <- 1/sqrt(Nsim)
Bart_se_pacf_rest <- sqrt((1+2*cumsum(test2_pacf$acf^2)[-Nlag])/Nsim)
Bart_se_pacf <- c(Bart_se_pacf_frist, Bart_se_pacf_rest)
ConfBands <- qnorm((1-0.05/2)) * Bart_se_pacf
lines(test2_pacf$lag, ConfBands, col="red") #add the calculated CI to the plot, just to check that it goes in the right spot
lines(test2_pacf$lag, ConfBands*-1, col="red")

A: Here is the extract from the p.16 and corresponding code from p. 17 of the book by Bernhard Pfaff which answers your question: 
The assumption of uncorrelatedness can be tested with the Ljung-Box Portmanteau test (see Ljung and Box [1978]). This test is implemented in the Box.test() function of the package stat.Except for the PACF, these tools are graphically returned by the function tsdiag(). Note in the following code y is the variable of an interest (e.g, GDP or unemployment rate). 
## ARMA ( 1 , 1 )
arma11 <− arima(y , order = c(1, 0, 1) )
tsdiag(arma11)
res11 <− residuals (arma11)
Box . test ( res11 , lag = 20, type = ”Ljung−Box ”)
tsdiag(arma11)

A: I have written an implementation that is identical to STATA. plot = TRUE gives a plot. Happy if someone maintaining a time-series package would like to include it. 
corrgrm <- function(x, lag.max = NULL, plot = FALSE) {
  if(plot) {
    oldpar <- par(mfrow = c(1, 2))
    on.exit(par(oldpar))
  }
  ac <- eval(substitute(acf(x, lag.max, plot = plot)))
  acf <- drop(ac[[1L]])[-1]
  pacf <- drop(eval(substitute(pacf(x, lag.max, plot = plot)))[[1L]])
  n <- ac$n.used
  # n * cumsum(acf^2) Box Pierce -> poor small sample properties
  qstat <- n * (n + 2) * cumsum(acf^2 / seq.int(n - 1, n - length(acf))) # Ljung-Box: The standard: Same as STATA !!
  lags <- seq_along(acf)
  res <- cbind(Lag = lags, AC = acf, PAC = pacf,
               Q = qstat, `Pr(>Q)` = pchisq(qstat, lags, lower.tail = FALSE))
  class(res) <- "corrgrm"
  res
}

print.corrgrm <- function(x, digits = 3, ...) {
  xx <- cbind(format(x[, 1L, drop = FALSE], drop0trailing = TRUE),
              format(round(x[, -1L], digits)))
  print.default(`rownames<-`(xx, rep(" ", nrow(xx))), quote = FALSE, right = TRUE)
}

