Solved – 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?

Contents

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

Rate this post