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?
 A: The problem is solved once you can solve it for lag 1, because you can take $k$ such independent (or at least uncorrelated) time series $X_t^{(1)},$ $X_t^{(2)},$ through $X_t^{(k)}$ and interleave them to form a time series
$$X_t = X_1^{(1)}, X_1^{(2)}, \ldots, X_1^{(k)},\ X_2^{(1)}, \ldots, X_2^{(k)},\ X_3^{(1)}, \ldots$$
with zero autocorrelation at all but lags that are multiples of $k.$  Let's call this the "$k$-interleaved series."
Consider a portion of a time series process $X_t$, which we may index at times $t=1,2,\ldots,n.$  For the autocorrelations to be meaningful, we must assume the entire process is at least weakly second order stationary.  The question requires the covariance matrix of the random vector $\mathcal{X}=(X_1,X_2,\ldots,X_n)$ to have the form
$$\mathbb{P}_n(\rho) = \pmatrix{1 & \rho & \color{gray}0 & \color{gray}0 & \cdots & \color{gray}0 &\color{gray}0\\
\rho & 1 & \rho & \color{gray}0 & \cdots & \color{gray}0& \color{gray}0\\
\color{gray}0 & \rho & 1 & \rho & \cdots & \color{gray}0& \color{gray}0\\
\vdots & \vdots & \vdots & \ddots & \cdots & \vdots& \vdots\\
\color{gray}0 & \color{gray}0 & \color{gray}0 & \color{gray}0 & \cdots  &\rho & \color{gray}0\\
\color{gray}0 & \color{gray}0 & \color{gray}0 & \color{gray}0 & \cdots  &1& \rho\\
\color{gray}0 & \color{gray}0 & \color{gray}0 & \color{gray}0 & \cdots & \rho & 1
}$$
where $\rho$ is the common lag-1 autocovariance.
Regardless what value $\rho$ might have, this matrix has the eigenvectors
$$q_{j;n} = \left(\sin\left(\frac{2\pi i j}{n+1}\right)\right)_{i=1,2,\ldots,n}$$
with corresponding eigenvalues $\lambda_{j;n}(\rho)$ (which may easily be computed from the eigenvalue equation: $\lambda_{j;n}$ is the common ratio of the components of $\mathbb{P}_n(\rho)q_{j;n}$ and the components of $q_{j;n}$).  Provided $|\rho|\le 1/2,$ all those eigenvalues will be non-negative, which implies there exists a random vector $\mathcal X$ with $\mathbb{P}_n(\rho)$ as its covariance matrix.  The orthogonal matrix $\mathbb{Q}_n = (q_{ij}) = (q_1;q_2;\cdots;q_n)$ whose columns comprise the normalized eigenvectors diagonalizes $\mathbb{P}_n(\rho).$

Consequently, when $\mathcal{Y}_n(\rho) = (Y_{j;n}(\rho))$ is any
vector of $n$ independent zero-mean variables with variances
$$\operatorname{Var}(Y_{j;n}) = \lambda_{j;n}(\rho),$$ the variable
$$\mathcal{X}_n(\rho) = \mathbb{Q}_n \mathcal{Y}_n(\rho)$$
has all the properties desired of $X_t.$


This demonstration leads to a simple and fairly efficient algorithm to generate realizations of the process $X_t.$  (It requires $O(n^2)$ computation, but likely could be sped up to $O(n\log(n))$ using the FFT.)  On this workstation it requires less than 0.02 seconds to generate one realization with $n=2000.$
Here are the empirical acfs of eight simulated series with $n=2000$ and lag-1 coefficient $\rho=1/2:$

All the lags but the first differ only randomly from zero, while the lag-1 coefficient differs only randomly from $\rho.$
When these $8$ series are interleaved, they form one series of length $8\times 2000=16000$ and, by design, the only nonzero autocorrelation coefficient occurs at lag $8,$ where it equals $\rho:$

What would such a series look like?  Here are the first $128$ values of the interleaved series from the simulation:

There is a hint of a seasonal cycle of period 8.  However, as we know, this series has no seasonal fluctuations.  For instance, the stl function in R finds only a tiny seasonal component of variance $0.0017$ compared to residuals of variance $0.82,$ almost 500 times larger.
This is the R code used to produce the simulated data and graphics.
n <- 2000     # Time series length (2 or longer)
rho <- 0.5    # Lag-1 autocorrelation (safe between -0.5 and +0.5).
n.sim <- 8    # Number of realizations to create
rf <- rnorm   # Generates zero-mean random values
system.time({
  #
  # Compute eigenvector components.
  #
  y <- sin(seq(0, 2*pi, length.out=2*n+3))[-(2*n+3)]
  y <- y / sqrt(sum(y^2)/2) # Normalized
  #
  # Compute eigenvalues.
  #
  f <- function(i, y, rho) {
    n1 <- length(y) / 2
    j <- 0:2 * i %% (2*n1) + 1
    sum(y[j] * c(rho, 1, rho)) / y[j[2]]
  }
  lambda <- sapply(1:n, f, y=y, rho=rho)
  if (min(lambda) < 0) warning("Correlation ", rho, " for n = ", n, " is impossible.")
  #
  # Generate random values.
  #
  e <- matrix(rf(n * n.sim), n) * sqrt(lambda) # In columns
  #
  # Take linear combinations of them.
  #
  x <- t(sapply(seq_len(n), function(i) {
    j <- (1:n * i) %% (2*(n+1)) + 1
    colSums(e * y[j])
  }))
})
#
# The interleaved ACF.
#
acf(c(t(x)), lwd=2, main=bquote(paste("ACF of the ", .(n.sim), "-interleaved series")))
#
# Plot the ACFs.
#
par(mfrow=c(2,4))
invisible(apply(x, 2, acf, lag.max=10, main="", bty="n", lwd=2))
par(mfrow=c(1,1))
#
# Plot part of one simulated series.
#
plot(ts(c(t(x))[1:128]), ylab="Simulated Values", 
     main=bquote(paste("Prefix of the ", .(n.sim), "-interleaved series for ", rho==.(rho))))

A: Edit: It looks like I misinterpreted this question.  The question appears to be asking for a method to generate a random vector with non-zero auto-correlation at one lag.  My answer give a method to generate a random vector with non-zero auto-regression at only one lag, which is a different thing (and gives non-zero auto-correlation at multiple lags).  I am going to leave the answer here because I think it may still be useful to some users interested in the content anyway.

Use the rGARMA function in the ts.extend package
You can generate random vectors from any stationary Gaussian ARMA model using the ts.extend package.  This package generates random vectors directly form the multivariate normal distribution using the computed autocorrelation matrix for the random vector, so it gives random vectors from the exact distribution and does not require "burn-in" iterations.  Here is an example from an $\text{AR}(2)$ model with non-zero autoregression only at lag $\ell = 2$.
#Load the package
library(ts.extend)

#Set parameters
AR    <- c(0, 0.4)
m     <- 100

#Generate n = 12 random vectors from this model
set.seed(1)
SERIES <- rGARMA(n = 12, m = m, ar = AR)

#Plot the series using ggplot2 graphics
library(ggplot2)
plot(SERIES)


