Here is such an example. I will describe the process in terms of its sample paths.
It's a simple process, it only has four sample paths: $\omega_1$, $\omega_2$, $\omega_3$, $\omega_4$. This means that the only possible outcomes are:
$$ ...,\ X(-1) = \omega_i(-1),\ X(0) = \omega_i(0),\ X(1) = \omega_i(1),\ ... $$
One outcome for every $i \in \{1,2,3,4\}$.
The sample paths are defined as follows, for $t\in \mathbb{Z}$:
$$\omega_1(t):= \begin{cases}
1 & \text{if $t$ is a multiple of $4$,}\\
0 & \text{else.}\\
\end{cases} $$
$$\omega_2(t):= \begin{cases}
0 & \text{if $t$ is a multiple of $4$,}\\
1 & \text{else.}\\
\end{cases} $$
$$\omega_3(t):= \omega_1(t-2) $$
$$\omega_4(t):= \omega_2(t-2) $$
So they are:
$$ ...,1,0,0,0,1,0,0,0,1,0,0,0,1,0,0,0,... $$
$$ ...,0,1,1,1,0,1,1,1,0,1,1,1,0,1,1,1,... $$
$$ ...,0,0,1,0,0,0,1,0,0,0,1,0,0,0,1,0,... $$
$$ ...,1,1,0,1,1,1,0,1,1,1,0,1,1,1,0,1,... $$
Here is a plot of the four possible sample paths:
Each sample path is defined to have probability $1/4$. This concludes the definition of the process.
I will now show that it satisfies the conditions required by the OP.
For $t\in\{1,2,3,4\}$ one can see that there are exactly $2$ sample paths such that $X(t)=1$, and exactly $2$ sample paths such that $X(t)=0$. Hence for $t\in\{1,2,3,4\}$:
$$ P[X(t)=1] = P[X(t)=0] = 1/2 .$$
Since the process is periodic this is true for every $t\in \mathbb{Z}$. This means that the distribution of $X(t)$ does not depend on $t$ and therefore the process is 1st order stationary.
Likewise, looking at consecutive pairs, for $t\in\{1,2,3,4\}$ there is exactly $1$ sample path with $X(t)=1$ and $X(t+1)=0$, exactly $1$ sample path with $X(t)=0$ and $X(t+1)=1$, exactly $1$ sample path with $X(t)=0$ and $X(t+1)=0$, and exactly $1$ sample path with $X(t)=1$ and $X(t+1)=1$. So for $t\in\{1,2,3,4\}$:
$$ P[X(t)=1,X(t+1)=0] = P[X(t)=0,X(t+1)=1] =P[X(t)=0,X(t+1)=0] =P[X(t)=1,X(t+1)=1] = 1/4 .$$
Again by periodicity, this is true for every $t\in \mathbb{Z}$ and the process is 2nd order stationary.
However, for instance: $$P[X(0)=0,X(1)=0,X(2)=0] = 0,$$ while $$P[X(1)=0,X(2)=0,X(3)=0] = 1/4.$$
Therefore the process is not 3rd order stationary
Second order stationarity is weaker than strict stationarity. Second order stationarity requires that first and second order moments (mean, variance and covariances) are constant throughout time and, hence, do not depend on the time at which the process is observed. In particular, as you say, the covariance depends only on the lag order, $k$, but not on the time at which it is measured, $Cov(x_t, x_{t-k}) = Cov(x_{t+h}, x_{t+h-k})$ for all $t$.
In a strict stationarity process, the moments of all orders remain constant throughout time, i.e., as you say, the joint distribution of $X_{t1},X_{t2},...,X_{tm}$ is the same as the joint distribution of $X_{t1+k}+X_{t2+k}+...+X_{tm+k}$ for all $t1,t2,...,tm$ and $k$.
Therefore, strict stationarity involves second order stationarity but the converse is not true.
Edit (edited as answer to @whuber's comment)
The previous statement is the general understanding of weak and strong stationarity. Although the idea that stationarity in the weak sense does not imply stationarity in a stronger sense may agree with intuition, it may not be so straightforward to proof, as pointed out by whuber in the comment below. It can be helpful to illustrate the idea as suggested in that comment.
How could we define a process that is second-order stationary (mean, variance and covariance constant throughout time) but it is not stationary in strict sense (moments of higher order depend on time)?
As suggested by @whuber (if I understood correctly) we can concatenate batches of observations coming from different distributions. We just need to be careful that those distributions have the same mean and variance (at this point let's consider that they are sampled independently of each other). On one hand, We can for example generate observations from the Student's $t$-distribution with $5$ degrees of freedom. The mean is zero and the variance is $5/(5-2)=5/3$. On another hand, we can take the Gaussian distribution with zero mean and variance $5/3$.
Both distributions share the same mean (zero) and variance ($5/3$). Thus, the concatenation of random values from these distribution will be, at least, second-order stationary. However, the kurtosis at those points governed by the Gaussian distribution will be $3$, while at those time points where the data come from the Student's $t$-distribution it will be $3+6/(5-4)=9$. Therefore, the data generated in this way are not stationary in strict sense because moments of fourth order are not constant.
The covariances are also constant and equal to zero, since we considered independent observations. This may seem trivial, so we can create some dependence among observations according the following autoregressive model.
$$
y_t = \phi y_{t-1} + \epsilon_t \,, \quad |\phi| < 1 \,, \quad t = 1,2,...,120
$$
with
\begin{eqnarray}
\epsilon_t \sim \left\{
\begin{array}{ll}
N(0, \sigma^2=5/3) \quad & \hbox{if} \; t \in [0,20], [41,60], [81,100] \\
t_5 \quad & \hbox{if} \; t \in [21,40], [61,80], [101,120] \,.
\end{array}
\right.
\end{eqnarray}
$|\phi| < 1$ ensures that second-order stationarity is satisfied.
We can simulate some of these series in the R software and check whether the sample mean, variance, first order covariance and kurtosis remain constant across batches of $20$ observations (the code below uses $\phi=0.8$ and sample size $n=240$, the Figure displays one of the simulated series):
# this function is required below
kurtosis <- function(x)
{
n <- length(x)
m1 <- sum(x)/n
m2 <- sum((x - m1)^2)/n
m3 <- sum((x - m1)^3)/n
m4 <- sum((x - m1)^4)/n
b1 <- (m3/m2^(3/2))^2
(m4/m2^2)
}
# begin simulation
set.seed(123)
n <- 240
Mmeans <- Mvars <- Mcovs <- Mkurts <- matrix(nrow = 1000, ncol = n/20)
for (i in seq(nrow(Mmeans)))
{
eps1 <- rnorm(n = n/2, sd = sqrt(5/3))
eps2 <- rt(n = n/2, df = 5)
eps <- c(eps1[1:20], eps2[1:20], eps1[21:40], eps2[21:40], eps1[41:60], eps2[41:60],
eps1[61:80], eps2[61:80], eps1[81:100], eps2[81:100], eps1[101:120], eps2[101:120])
y <- arima.sim(n = n, model = list(order = c(1,0,0), ar = 0.8), innov = eps)
ly <- split(y, gl(n/20, 20))
Mmeans[i,] <- unlist(lapply(ly, mean))
Mvars[i,] <- unlist(lapply(ly, var))
Mcovs[i,] <- unlist(lapply(ly, function(x)
acf(x, lag.max = 1, type = "cov", plot = FALSE)$acf[2,,1]))
Mkurts[i,] <- unlist(lapply(ly, kurtosis))
}
The results are not what I expected:
round(colMeans(Mmeans), 4)
# [1] 0.0549 -0.0102 -0.0077 -0.0624 -0.0355 -0.0120 0.0191 0.0094 -0.0384
# [10] 0.0390 -0.0056 -0.0236
round(colMeans(Mvars), 4)
# [1] 3.0430 3.0769 3.1963 3.1102 3.1551 3.2853 3.1344 3.2351 3.2053 3.1714
# [11] 3.1115 3.2148
round(colMeans(Mcovs), 4)
# [1] 1.8417 1.8675 1.9571 1.8940 1.9175 2.0123 1.8905 1.9863 1.9653 1.9313
# [11] 1.8820 1.9491
round(colMeans(Mkurts), 4)
# [1] 2.4603 2.5800 2.4576 2.5927 2.5048 2.6269 2.5251 2.5340 2.4762 2.5731
# [11] 2.5001 2.6279
The mean, variance and covariance are relatively constant across batches as expected for a second-order stationary process. However, the kurtosis remains relatively constant as well. We could have expected higher values of the kurtosis at those batches related to draws from the Student's $t$-distribution. Maybe $20$ observations is not enough to capture changes in kurtosis. If we didn't know the data generating process of these series and we looked at rolling statistics, we would probably conclude that the series is stationary at least up to order fourth. Either I didn't take the right example or some features of the series get masked for this sample size.
Best Answer
It might help to think about it in terms of copulae. $X$ and $Y$ can be converted to grades and joined together by a copula. This copula could potentially be constant or a time varying DCC copula. As a result, the joint distribution could be potentially multivariate constant in the first and second moments or time-varying. In that sense, they would be univariate stationary, but not necessarily multivariate stationary unless you make some assumptions about the dependence structure.