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?

**Contents**hide

#### Best Answer

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") `

### Similar Posts:

- Solved – How many lags for ADF test based on ACF, PACF
- Solved – say the residuals have an ARCH effect in this plot
- Solved – How are the p-values in the glmer (lme4) output calculated
- Solved – what does it mean if the partial auto correlation function of a time series have a value >+1 or <- 1
- Solved – what does it mean if the partial auto correlation function of a time series have a value >+1 or <- 1