Let me comment on points 4) and 5). The problem with infinities in QFT or traditional equilibrium statistical field theory
is related to the one addressed by Martin's theory but there are some differences. For concreteness let me talk about the $\phi^4$ model
only. Mathematically, the problem it poses is to make sense of the probability measure
$$
\frac{1}{\mathcal{Z}}\exp\left(
-\int_{\mathbb{R}^d}\{
\frac{1}{2} (\nabla\phi)^2(x)+\mu \phi(x)^2+g\phi(x)^{4}
\} d^dx
\right)\ D\phi
$$
on the "space of all functions" $\phi:\mathbb{R}^d\rightarrow\mathbb{R}$. This is the kind of heuristic formulas
one finds in physics QFT textbooks. The symbol $D\phi$ stands for Lebesgue measure on this space of functions
and $\mathcal{Z}$ is a normalization constant so the full space has measure one as befits a probability measure.
Now let's turn this into a well posed mathematical question.

First remove the $\phi^2$ and $\phi^4$ terms, i.e., consider the case $\mu=g=0$. Then this measure $\mu_{C_{-\infty}}$
makes perfect sense. It is the centered Gaussian measure on the space of temperate distributions $S'(\mathbb{R}^d)$
and with covariance $C_{-\infty}$ given by
$$
C_{-\infty}(f,g)=\frac{1}{(2\pi)^{d}}\int_{\mathbb{R}^d}\frac{\overline{\widehat{f}(\xi)}
\widehat{g}(\xi)}{|\xi|^{2}} d^d\xi
$$
for all test functions $f$ and $g$ in $S(\mathbb{R}^d)$. Using this first rigorous step, one can reformulate 
the problem as that of making sense of 
$$
\frac{1}{\mathcal{Z}}\exp\left(
-\int_{\mathbb{R}^d}\{
\alpha (\nabla\phi)^2(x)+\mu \phi(x)^2+g\phi(x)^{4}
\} d^dx
\right)\ d\mu_{C_{-\infty}}(\phi)
$$
with a new normalization constant $\mathcal{Z}$ that I will still keep calling $\mathcal{Z}$.
I also introduced the ``wave function renormalization coupling constant'' $\alpha$ for more generality.
We made a bit of progress (we avoided the problematic Lebesgue measure $D\phi$),
but this still does not make mathematical sense because $\mu_{C_{-\infty}}$ is supported on nasty Schwartz distributions
and pointwise powers like $\phi^2$ and $\phi^4$ are ill-defined, as in Martin's SPDE setting. This is the source
of the UV (ultraviolet) infinities. There are also IR (infrared) problems due to the integration inside the exponential
being over $\mathbb{R}^d$ instead of a compact set. To address these issues, we need what the French call *troncature et
régularisation*. Let $\rho_{\rm UV}$ be a mollifier, i.e., a compactly supported $C^{\infty}$ function
$\mathbb{R}^d\rightarrow\mathbb{R}$ with $\int \rho_{\rm UV}=1$.
Let $\rho_{\rm IR}$ be a cut-off function, i.e., a nonnegative compactly supported $C^{\infty}$ function
$\mathbb{R}^d\rightarrow\mathbb{R}$ which is equal to 1 in a neighborhood of the origin. To slice Fourier momenta
into (Littlewood-Paley) shells we introduce an integer $L>1$, not necessarily equal to 2 as is customary in
harmonic analysis. For $r,s\in\mathbb{Z}$, define the rescaled functions $\rho_{{\rm UV},r}(x)=L^{-dr}\rho_{\rm UV}(L^{-r}x)$
and $\rho_{{\rm IR},s}(x)=\rho_{\rm IR}(L^{-s}x)$, and consider the probability measure $\nu_{r,s}$ given by
$$
\frac{1}{\mathcal{Z}}\exp\left(
-\int_{\mathbb{R}^d}\rho_{{\rm IR},s}(x)\{
\alpha (\nabla\phi)^2(x)+\mu \phi(x)^2+g\phi(x)^{4}
\} d^dx
\right)\ d\mu_{C_{r}}(\phi)
$$
where $\mu_{C_r}$, or regularized Gaussian measure, is the direct image of $\mu_{C_{-\infty}}$
by the convolution map $\phi\mapsto \rho_{{\rm UV},r}\ast\phi$.
In other words, $\mu_{C_r}$ is the centered Gaussian measure with covariance
$$
C_{r}(f,g)=\frac{1}{(2\pi)^{d}}\int_{\mathbb{R}^d}\frac{|\widehat{\rho_{{\rm UV},r}}(\xi)|^2\ \overline{\widehat{f}(\xi)}
\widehat{g}(\xi)}{|\xi|^{2}} d^d\xi\ .
$$
A good metaphor would be to say that your orginal flat-screen TV was too smart. The linear size of the
screen was $L^s=\infty$ and that of a pixel was $L^r=0$. Instead one should make $r$ and $s$ finite so that $\nu_{r,s}$
is mathematically well defined, and then study the limit where $r\rightarrow-\infty$ and $s\rightarrow\infty$
in the sense of
weak convergence of probability measures on the topological space $S'(\mathbb{R}^d)$.
Now renormalization theory in physics tells us that unless we allow the couplings $(\alpha,\mu,g)$ to depend on the
UV cut-off scale $r$, the following is more likely to happen: 1) we don't converge (e.g., loss of tightness), 2)
we converge to something utterly uninteresting like the atomic measure on the singleton $\{\phi=0\}$, 3) we
converge to something less trivial but still uninteresting, namely, a Gaussian measure like the GFF $\mu_{C_{-\infty}}$
or white noise, or massive free fields interpolating between the two.
Therefore the weak limit we need to study depends on the choice of bare ansatz $(\alpha_r,\mu_r,g_r)_{r\in\mathbb{Z}}$
(or rather the germ of this sequence at $r=-\infty$). Finally the well posed mathematical question I promised,
regarding trying to mense of the original $\phi^4$ functional integral is the following.

**Problem:** Find an explicit parametrization of all weak limits (of probability measures on $S'(\mathbb{R}^d)$)
given by $\lim_{r\rightarrow-\infty}\lim_{s\rightarrow\infty}\nu_{r,s}$ for all possible choices of bare ansatz
$(\alpha_r,\mu_r,g_r)_{r\in\mathbb{Z}}$.

What renormalization theory in physics also tells us is that although it seems one has a hugely infinite-dimensional
amount of freedom in choosing the sequence $(\alpha_r,\mu_r,g_r)_{r\in\mathbb{Z}}$, the set $\mathscr{T}$
of weak limit points
is a finite dimensional variety. For $d=3$, one expects three parameters or "renormalized coupling constants"
$(\alpha_{\rm R},\mu_{\rm R},g_{\rm R})$ suffice.
One can even get rid of $\alpha_{\rm R}$ if one quotients by taking constant
multiples of the random field $\phi$. 

There are rigorous renormalization group techniques for constructing portions of $\mathscr{T}$. Kupiainen's work mentioned
by Ofer is an adaptation of these techniques to the time-dependent rough SPDE setting. 
The above is what I call the *non-anchored* Gibbsian way
of trying to construct elements $\nu\in\mathscr{T}$. There is a completely different approach which I call the
*anchored* stochastic quantization approach. There one need to make sense of the SPDE in Martin's answer which he
did, locally in time. Then one needs to understand this SPDE globally in time and construct an invariant
measure which gives $\nu\in\mathscr{T}$. This is also fraught with difficulties but there has been some progress
(see, e.g., <a href="https://arxiv.org/abs/1610.03415">this article</a> related to invariant measures).
The key difference between the two approches is that in the non-anchored setting one does not have a fixed probability space
to work with. In the second anchored situation one does since all fields are functionals of the driving noise. In the anchored setting,
$L^2$ estimates involving second moments only are enough to prove convergence in probability and thus in
law for the random fields of interest. In the non-anchored situation one needs to control all moments
(correlation functions) with uniform $n!$ bounds on these moments.

It is hard to say more as an MO answer, but you can see these preprints for further explanations:

 1. <a href="https://arxiv.org/abs/1604.05259">"A second-quantized Kolmogorov-Chentsov theorem"</a>
 2. <a href="https://arxiv.org/abs/1511.03180">"Towards three-dimensional conformal probability"</a>
 3. <a href="https://arxiv.org/abs/1311.4897"> "QFT, RG, and all that, for mathematicians, in eleven pages"</a>

The article in bullet point 1) provides an alternate way of defining pointwise products of random distributions
(like $\phi^2$ and $\phi^4$ above) using Wilson's operator product expansion (OPE). Regarding 5) in the OP's question,
I believe it would be of great interest to prove that the moments of the solution $\Phi$ constructed by Martin
satisfy the (dynamical version) of Wilson's OPE, and then compare $\Phi^3$ obtained by the theory of regularity structures
with the one constructed (from the OPE) in my 2nd quantized KC paper.