How can I generate a time series with autocorrelation at lags other than 1?

How can I generate a time series that has autocorrelation at a certain lag, but only that lag and nothing else?

If you don't require stationarity for your process then generating it is as simple as starting with any arbitrary random variable and then using the recursive equations of the desired process to generate the rest of the series. However, assuming that you want to generate a stationary process, this requires you to find the stationary distribution. By way of example, I will show you how to generate a stationary $$\text{AR}(2)$$ process with missing coefficient on its first-order lag.

Deriving the stationary distribution: The $$\text{AR}(2)$$ process of interest here is:

$$X_t = \phi X_{t-2} + \varepsilon_t \quad \quad \quad \quad \quad \varepsilon_t \sim \text{IID N}(0, \sigma).$$

(Note that I have only included an auto-regressive term on the last lag, as was specified in your question.) Stationarity of the process requires $$|\phi| < 1$$. With a bit of algebra, the Wold representation for the process can be written as:

$$X_t = \sum_{k=0}^\infty \psi_k \cdot \varepsilon_{t-k} \quad \quad \quad \quad \quad \psi_k = \begin{cases} \phi^{k/2} & & & \text{for even } k, \\[6pt] 0 & & & \text{for odd } k. \\ \end{cases}$$

Using the Wold representation we get the variance:

\begin{equation} \begin{aligned} \mathbb{V}(X_t) &= \sum_{k=0}^\infty \psi_k^2 \cdot \mathbb{V}(\varepsilon_{t-k}) \\[6pt] &= \sum_{k=0}^\infty \phi^k \cdot \mathbb{V}(\varepsilon_{t-2k}) \\[6pt] &= \sigma^2 \sum_{k=0}^\infty \phi^k \\[6pt] &= \frac{\sigma^2}{1-\phi}. \\[6pt] \end{aligned} \end{equation}

(We also have the covariance $$\mathbb{C}(X_t,X_{t-1}) =0$$.) Thus, the stationary distribution of the process is:

$$X_\infty \sim \text{N} \bigg( 0, \frac{\sigma^2}{1-\phi} \bigg).$$

Generating the process: Now that we have the stationary distribution of the process, we can generate a vector of values from the process by generating two consecutive values from the stationary distribution (noting that they are uncorrelated), and then generating the remaining values from the recursive equations. Here is some R code that generates a series of values from the above stationary model.

#Function to generate vector from the AR(2) process set out above
PROCESS <- function(T, phi = 0, sigma = 1) {
if (abs(phi) >= 1) { stop("Error: Process is non-stationary"); }
if (sigma < 0)     { stop("Error: Process standard error is negative"); }
VALUES      <- numeric(T);
ERROR       <- rnorm(T, mean = 0, sd = sigma);
VALUES[1:2] <- ERROR[1:2]/sqrt(1-phi);
if (T > 2) { for (t in 3:T) {
VALUES[t] <- phi*VALUES[t-2] + ERROR[t]  } }
VALUES; }

We can use this function to generate a vector of length $$T$$ from the above $$\text{AR}(2)$$ process, for any specification of $$\phi$$ and $$\sigma$$ (we restrict the former to require stationarity of the process).

#Set values for time-series
T     <- 100;
phi   <- 0.4;
sigma <- 1;

#Generate and plot the series
set.seed(1);
TTT <- PROCESS(T, phi, sigma);
plot(TTT, type = 'l', main = 'AR(2) time-series',
xlab = 'Time', ylab = 'Value'); 