Dealing with one-variable and smooth function $f$ on a real interval $I$ such that $D^m f\in\mathcal{C}^2$, we have by Taylor theorem centered at $a\in I$ $$ D^mf(y)= D^mf(a) + D^{m+1}f(a)(y-a) + \frac{1}{2} D^{m+2}f(a) (y-a)^2 + \underbrace{\phi_a(y)(y-a)^2}_{remainder}, $$ where $\phi_a$ is a continuous function with $\lim_{y\to a}\phi_a(y)=0$. Now integrating as follows $$ \int_a^x \int_a^{z_{m-1}}\cdots\int_a^{z_1} D^m f(y)\ dy\ dz_1...dz_{m-1} $$ gives $$ f(x) - f(a) - D^1f(a)(x-a) - ... -\frac{1}{(m-1)!}D^{m-1}f(a)(x-a)^{m-1}\quad .$$
I want to do this in $N$-dimensions. Given $f=f(x_1,...,x_N)$ and $a=(a_1,...,a_N)\in \Omega \subseteq\mathbb{R}^N$, assume $f\in\mathcal{C}^{m+2} ( \Omega)$, then Taylor theorem for $D^\alpha f:=\frac{\partial^{\vert \alpha\vert} f}{\partial x_1^{\alpha_1}...x_N^{\alpha_N}}$ reads now
$$ D^\alpha f (y) = D^\alpha f(a) + \nabla D^\alpha f(a) (y-a) + \frac{1}{2} (y-a)^t H_{D^\alpha f(a)} (y-a) + \sum_{\vert \alpha \vert = 2} \phi_a (y) (y-a)^2 \quad (*) , $$
where $\nabla g = (D^1_1 g , D^1_2 g,...,D^1_N g)$, $H_{D^\alpha f(a)} $ is the hessian matrix of $D^\alpha f$ evaluated at $a$, $\alpha = (\alpha_1,...,\alpha_N)\in\mathbb{N}^N$ is a mult-index with $|\alpha | = \alpha_1 + \cdots + \alpha_N$ and $(y-a)^\alpha = (y_1-a_1)^{\alpha_1}\cdots (y_N-a_N)^{\alpha_N} .$
Over what do I should integrate the expression $(*)$ in order to get the analogous result as in one dimension, $$ f(x) - f(a) - \sum_{\vert \alpha \vert \leq m-1} D^\alpha f(a) (y-a)^\alpha $$
Thanks in advance. Any comment is welcome! I asked this question twice at mathexchange, but i got no answers nor comments.