This problem falls to the class called structured-output learning. First, you need to define a function $F(\mathbf{y}; \mathbf{x}, \mathbf{w})$ that scores predictions. It is supposed to return large values for the correct $\mathbf{y}$ and those similar to it. It is parametrized by the vector $\mathbf{w}$; usually $F$ depends linearly on $\mathbf{w}$. At test time, the decision is thus made by maximization:
$\hat{\mathbf{y}} = \arg\!\max_{\bar{\mathbf{y}}} F(\bar{\mathbf{y}}; \mathbf{x}, \mathbf{w})$.
So $F$ should be chosen the way to perform this maximization efficiently. See below for an example.
At training time, you need to find the parameters $\mathbf{w}^*$ such that for each training example $\mathbf{x}_i$ the corresponding prediction $\hat{\mathbf{y}}_i$ should be close $\mathbf{y}_i$ in some sense. Similarly to regression, it is unlikely to match perfectly, so one needs to define a loss function.
One such method is provided by structured SVM. The methods tries to find parameters $\mathbf{w}$ to maximize the margin
$F(\mathbf{y}_i; \mathbf{x}_i, \mathbf{w}) - \max_{\bar{\mathbf{y}}: \bar{\mathbf{y}} \ne \mathbf{y}_i} F(\bar{\mathbf{y}}; \mathbf{x}_i, \mathbf{w})$
for each training example $(\mathbf{x}_i, \mathbf{y}_i)$. To account for similarity in $\mathbf{y}$'s, it is usually augmented with the empirical loss function $\Delta(\bar{\mathbf{y}}; \mathbf{y})$ that measures some sort of distance; for discrete vectors it may measure Hamming distance. Thus, for each training example we want to maximize the following:
$F(\mathbf{y}_i; \mathbf{x}_i, \mathbf{w}) - \max_{\bar{\mathbf{y}}} \{ F(\bar{\mathbf{y}}; \mathbf{x}_i, \mathbf{w}) + \Delta(\bar{\mathbf{y}}; \mathbf{y}_i) \}$.
The idea is to strive for the larger margin if $\bar{\mathbf{y}}$ is more different from $\mathbf{y}_i$. We can put it together and add $L_2$ regularization to get the SVM-style optimization problem:
$\textrm{minimize}_{\mathbf{w}, \boldsymbol{\xi}} \quad \frac{\|\mathbf{w}\|^2}{2} + C \sum_{i=0}^N \xi_i,$
$\textrm{s.t.} \quad F(\mathbf{y}_i; \mathbf{x}_i, \mathbf{w}) - \max_{\bar{\mathbf{y}}} \{ F(\bar{\mathbf{y}}; \mathbf{x}_i, \mathbf{w}) + \Delta(\bar{\mathbf{y}}; \mathbf{y}_i) \} \ge -\xi_i, \quad \forall i.$
It is not trivial to solve, but there are well-studied methods that perform it, including the cutting-plane algorithm and subgradient descent. For implementations, you can try SVMstruct and pystruct.
What is left is the design of the scoring function $F$. It really depends on the problem you solve, but often some kind of graphical model is useful. Given that your data are sequential, something like Hidden Markov Model or some its generalization may be useful. For example, you may consider the following “factorization”:
$F(\mathbf{y}; \mathbf{x}, \mathbf{w}) = \psi(y_{0}; \mathbf{x}, \mathbf{w}) + \sum_{j=1}^q \Big( \phi(y_{j-1}, y_{j}; \mathbf{w}) + \psi(y_{j}; \mathbf{x}, \mathbf{w}) \Big).$
It allows performing maximization over $\mathbf{y}$ in $\mathcal{O}(q)$ time using dynamic programming. The potential functions $\psi$ and $\phi$ should depend linearly on (different parts of) $\mathbf{w}$, but specific form depends on your problem (in particular, if the output variables are discrete or continuous). You also need to figure out the components of $\mathbf{x}$ the potential function depends on.
This form of scoring functions accounts for the “correlations” of neighbouring $y$'s, but it is possible to add higher-order connections at some computational cost (i.e. potentials like $\phi_4(y_{j-4}, y_{j}; \mathbf{w})$). SVMstruct library contains an example of training HMM.
Note that in case of such sequential models training can be done by maximizing likelihood, which leads to a different objective function, but in general case it is infeasible. See this paper for more details on probabilistic training of sequential models.
