Taylor's formula with remainder for vector-valued functions Let $f: \mathbb{R}^n \to  \mathbb{R}^n $. Does there exist a generalization of Taylor's Theorem with Lagrange Remainder for such a vector-valued function?

The short answer is "yes". The multivariable arguments to f require partial derivatives of f to determine the coefficients of the polynomial terms. The vector valued output leads to vectors of polynomials. Combining results in a mathematical structure analogous to the Taylor Polynomial, but called a "jet". The jet is used in differential geometry and a good introduction can be found on the Wikipedia Jet(mathematics) article.
The remainder term is often written as it is for the one variable case. But unpacking the generalization of the notation can be tricky.  Essentially, you have remainders in each coordinate of the vector output. Those remainders can be written as $$ f_i^{(k+1)}(\xi_i) {(x-x_0)^{\otimes(k+1)} \over (k+1)!}$$ for some $\xi_i$ in the neighborhood $U$ of $x_0$ you consider.  This formula looks very similar to the one dimensional case, but note that the powers of $(x-x_0)$ have been generalized -- as have the derivatives of $f$. For example, in 2D when $k=1$, you have remainder terms with $(x-x_0)^2$, $(x-x_0)(y-y_0)$, and $(y-y_0)^2$, and you have all possible $(k+1)$-order partial derivatives of $f$. Finally, $\xi_i$ can differ for each dimension $i$ of the output. As far as I know, there is (generally) no single point $\xi$ for which the remainder can be evaluated, but I don't have a counter example. Certainly, you can replace $f_i^{(k+1)}(\xi_i)$ with $$M=\max_U |f_i^{(k+1)}|$$ and get a bound on the Remainder.
Notation for the jet of order $k$ about $x_0$ for the function $f$ seems to be $(J^k_{x_0}f)$ for the polynomial part. Thus: 
$$f(x) = (J^k_{x_0}f)(x) + {R_{k+1}(x) \over (k+1)!} (x-x_0)^{\otimes(k+1)}$$
