How to contract a 4th rank Levi Civita tensor? I recently started tensors and my textbook mentions the following identity:
\begin{equation}
\varepsilon_{ijkl}\varepsilon_{mnpl} = \delta_{im}\delta_{jn}\delta_{kp}-\delta_{im}\delta_{jp}\delta_{kn}+\delta_{in}\delta_{jp}\delta_{km}\\ \hspace{1.7cm}-\delta_{in}\delta_{jm}\delta_{kp}+\delta_{ip}\delta_{jm}\delta_{kn}-\delta_{ip}\delta_{jn}\delta_{km}
\end{equation}
How do I prove this? For the 3 dimensional analog, i.e.,
\begin{equation}
\varepsilon_{ijk}\varepsilon_{mnk} = \delta_{im}\delta_{jn} - \delta_{in}\delta_{jm}
\end{equation}
I verified the contraction identity by writing down each case for all indices. But that is not a proof, just a verification, and I think verifying in a similar fashion for 4th dimension is going to be cumbersome.
Is there any general prescription to prove such contractions, possibly for any $n$-th  order, greater than 4? Or at least for 4th order if not $n$?
 A: 
... and I think verifying in a similar fasion for 4th dimension is going to be cumbersome.

You are right.
After three dimensions, the regular Kronecker delta $\delta^i_j$ simply stops being useful on its own. We need to transition to the Generalized Kronecker Delta,
$$\delta^{i_1\dots i_p}_{j_1\dots j_p}=\begin{cases}0 & \text{if any }i\text{s or }j\text{s are repeated} \\ +1 & \text{if }j_1\dots j_p\text{ is an even permutation of }i_1\dots i_p \\ -1 & \text{if }j_1\dots j_p\text{ is an odd permutation of }i_1\dots i_p \end{cases}$$
It may be written as a determinant,
$$\delta^{i_1\dots i_p}_{j_1\dots j_p}=\det\begin{bmatrix}\delta^{i_1}_{j_1}&\cdots &\delta^{i_1}_{j_p} \\ \vdots&\ddots & \vdots \\ \delta^{i_p}_{j_1}&\dots&\delta^{i_p}_{j_p}\end{bmatrix}$$
Higher-order versions of the Levi-Civita symbol are generally defined using the generalized Kronecker delta:
$$\epsilon_{i_1\dots i_n}\equiv \delta^{1\dots n}_{i_1\dots i_n}$$
It can be shown that
$$\epsilon^{i_1\dots i_s~k_1\dots k_n}\epsilon_{j_1\dots j_s~k_1\dots k_n}=n!~\delta^{i_1\dots i_s}_{j_1\dots j_s}$$
Which in our case gives us
$$\epsilon^{ijkl}\epsilon_{mnpl}=\delta^{ijk}_{mnp}  =\det\begin{bmatrix}\delta^i_m &\delta^i_n &\delta^i_p \\ \delta^j_m& \delta^j_n& \delta^j_p \\\delta^k_m &\delta^k_n &\delta^k_p \end{bmatrix}$$
(Thankfully there will be 6 terms, as expected!)
Now you can expand this determinant to your hearts content.
