I think lag selection is usually done with shrinkage priors? You specify a large number of lags, then use something like a horseshoe prior to strongly bias the coefficients towards zero. See for example here.. Example implementation in PyMC here. I also do a modified version here, which has been suggested to be more compatible with NUTS. I’m not sure mine is implemented correctly.
Why are you thinking about correlations between the lags and the forecasts? What would that tell you? If the data generating process is stationary, the auto-covariance matrix should not be a function of time, so the correlations between time t and t-h for any t and a fixed h is the same. Thus the auto-correlations between the forecasts and any lags ought to be the same as the same auto-correlations as the in-sample data I think.