To formalize and generalize dlnB's +1 answer a little (based on Hamilton's textbook):
Cointegration implies that the deviations from the equilibrium are $I(0)$. Hence, some mechanism must bring back deviations back to the long-run relationship. This idea is formalized using error correction models (ECM).
Assume the following $VAR(p)$
\begin{equation}\tag{1}\label{1}
y_t = \alpha + \Phi_1{y_{t-1}} + \ldots + \Phi_p{y_{t - p}} + \epsilon_t
\end{equation}
Using the lag operators we can write this as
$$(I - \Phi_1{L^1} - \ldots - \Phi_pL^p) \cdot y_t = \Phi(L) y_t = \alpha + \epsilon_t$$
Define
\begin{equation}\tag{2}\label{2}
\rho \equiv \Phi_1 + \Phi_2 + \ldots + \Phi_p
\end{equation}
and
\begin{equation}\tag{3}\label{3}
\zeta_s \equiv - [\Phi_{s + 1} + \Phi_{s + 2} + \ldots + \Phi_p]
\end{equation}
Rewrite $I - \Phi_1{L^1} - \ldots - \Phi_pL^p$ by adding and immediately subtracting the coefficient matrices of order $j+1$ to $p$ on the lag operator of order $j$. We get
$$
\begin{gathered}
I - [(\Phi_1 + \Phi_2 + \ldots + \Phi_p) - (\Phi_2 + \Phi_3 + \ldots + \Phi_p)]L \hfill \\
- [(\Phi_2 + \ldots + \Phi_p) - (\Phi_3 + \ldots + \Phi_p)]L^2 \hfill \\
- [{\Phi_{p-1}} + \Phi_p - \Phi_p]L^{p-1} - \Phi_pL^p \hfill \\
\end{gathered}
$$
Using \eqref{2} and \eqref{3} yields
$$I - (\rho + \zeta_1)L - (\zeta_2 - \zeta_1)L^2 - \ldots - (\zeta_{p-1} - {\zeta_{p-2}})L^{p-1} - ( - \zeta_{p-1})L^p_ \cdot $$
Solving the terms in brackets gives
\begin{equation}\tag{4}\label{4}
I - \rho L - \zeta_1L - \zeta_2L^2+\zeta_1L^2 - \ldots - \zeta_{p-1}L^{p-1} + {\zeta_{p-2}}L^{p-1} - ( - \zeta_{p-1})L^p
\end{equation}
The $\zeta_i$-matrices appear both before the $i$th lag operator and, with reverse sign, before the $i+1$th lag operator. We can hence rewrite \eqref{4} as
$$I - \rho L - (\zeta_1L + \zeta_2L^2 + \ldots + \zeta_{p-1}L^{p-1})(1-L)$$
Hence, we have rewritten \eqref{1} as
$$\left[ I - \rho L - \left( \zeta_1L + \zeta_2L^2 + \ldots + \zeta_{p-1}L^{p-1} \right)(1-L) \right]y_t = \alpha + \epsilon_t$$
Multiplying out the square brackets, using $\Delta=1-L$, applying the lag operators and rearranging yields
$$y_t = \alpha + \rho y_{t-1} + \zeta_1\Delta y_{t-1} + \zeta_2\Delta y_{t-2} + \ldots + \zeta_{p-1}\Delta y_{t-p+1} + \epsilon_t$$
Subtract $y_{t-1}$ from either side to get
\begin{equation}\tag{5}\label{5}
\Delta y_t = \alpha - (I - \rho )y_{t-1} + \zeta_1\Delta y_{t-1} + \zeta_2\Delta y_{t-2} + \ldots + \zeta_{p-1}\Delta y_{t-p+1} + \epsilon_t
\end{equation}
Note $I-\rho=\Phi(1)$. This matrix is of reduced rank, say $h$. To see this, note the $I(1)$ assumption on $y_t$ implies that $\Phi(L)$ has a unit root, i.e.
$$
|I - \Phi_1{1^1} - \ldots - \Phi_p1^p|=|I - \Phi_1 - \ldots - \Phi_p|=0
$$
(note the determinant is 0 as matrix does not have full rank). That is, $\Phi(1)$ can be decomposed into two $(n\times h)$ matrices $B$ and $A$ such that
$$\Phi(1) = BA'$$
The $h$ rows of $A'$ are the cointegrating relationships. Linear combinations of cointegrating vectors and variables $A'y_{t-1}=e_{t-1}$ are stationary. We can thus rewrite \eqref{5} as
\begin{equation}\tag{6}\label{6}
\Delta y_t = \alpha + \zeta_1\Delta y_{t-1} + \zeta_2\Delta y_{t-2} + \ldots + \zeta_{p-1}\Delta y_{t-p+1}-Be_{t-1}+\epsilon_t
\end{equation}
Estimating the VAR in first differences implies omitting $Be_{t-1}$, which is relevant under cointegration.
