To understand whether your steps are correct you should compare to the mathematical formulation.
Step 1 is supposedly solving the rendering equation for some collection of initial rays
$$L_o(x_0, \omega_o) = L_e(x, \omega_o) + \int_{\mathcal{H}^2(x_0, N_{x_0})} f(\omega_o, x_0, \omega_i)L_i(x_0,\omega_i)(N_{x_0}\cdot \omega_i)\,d\omega_i.$$
Here $L_o(x_0,\omega_o)$ is the outgoing radiance function, i.e. the radiance leaving some surface point $x_0$ in direction $\omega_o$. $L_e(x_0,\omega_o)$ is the emitted radiance function, i.e. the radiance emitted from point $x_0$ in direction $\omega_o$. Next we have an integral over the hemisphere at point $x_0$ whose up direction is aligned with the normal $N_{x_0}$ at $x_0$. You can formalize this set as $\mathcal{H}^2(x_0,N_{x_0}) = \{z\in\mathbb{R}^3\,:\, \|z-x_0\|=1, \, (z-x_0)\cdot N_{x_0}\geq 0\}$, in practice one parametrizes it using spherical coordinates. $f$ is the BRDF. $L_i$ is the incident radiance function, i.e. $L_i(x_0,\omega_i) = L_o(r(x_0, \omega_i), -\omega_i)$, where $r$ is the ray-tracing function. Finally $(N_{x_0}\cdot \omega_i) = \cos\theta_i$ is the Lambertian cosine term that accounts for the angle between the ray and the surface.
For simplicity the above can be rewritten using linear operator notation as $L = L_e + TL$, from where one may get the Neumann expansion $L = \sum_{k=0}^{\infty}T^kL_e$, which is essentially a sum of increasingly dimensional integrals. You typically apply Monte Carlo to get an approximate solution of the latter.
The question is also how you choose the rays. Since you mentioned that you have an aperture, a film, and an exposure time, let us model the aperture as some set $\mathcal{A}\subset\mathbb{R}^3$, the film as $\mathcal{F}\subset\mathbb{R}^3$, and the exposure time as an interval $[t_0,t_1]$. Both of $\mathcal{F}$ and $\mathcal{A}$ are two-dimensional manifolds, and usually $\mathcal{A}$ is chosen to be a disk some distance in front of the film $\mathcal{F}$, and $\mathcal{F}$ is chosen to be a rectangle (in order to match the screen). The response at a given pixel (i,j) can be modeled as a weighted integration over the arriving radiance on the film, with some sensitivity function $w_{ij}$:
$$I_{ij} = \int_{t_0}^{t_1}\int_{\mathcal{F}}\int_{\mathcal{A}_{\mathcal{H}^2(x_{-1}, N_{x_{-1}})}}w_{ij}(x_{-1}, -\omega_o)L_i(x_{-1}, -\omega_o))(N_{x_{-1}}\cdot (-\omega_o))\, d(-\omega_o)dx_{-1}dt.$$
Here $\mathcal{A}_{\mathcal{H}^2(x_{-1}, N_{x_{-1}})}$ is the solid angle set of the aperture as seen from pixel $x_{-1}$. Here I have assumed that nothing depends on time (i.e. brdf, scene geometry, lights, aperture, film, sensitivity, all stay the same). With this assumption you can pop out the time integral as a $(t_1-t_0)$ multiplicative term as you describe in Step 5.
Note that physically the aperture is typically in front of the sensor, so I have modeled that instead of what you did. One could of course set the aperture behind the sensor in CG for convenience, and then you would get the solid angle of the film pixel's set as you do. I will however stick to the physical formulation of the aperture being in front of the sensor. We also need to get rid of that cosine dot product and want to integrate only over the pixel and not the whole film, which can be done by setting the sensitivity $w_{ij}(x_{-1},-\omega_o) = \frac{\boldsymbol{1}_{\mathcal{P}_{ij}}(x_{-1})}{N_{x_{-1}}\cdot (-\omega_o)}$. Here $\mathcal{P}_{ij}\subset \mathcal{F}$ is the set of points of the pixel on the film, and $\boldsymbol{1}_{\mathcal{X}}$ is the indicator function (it is one on the set, and zero outside of it). You want to average the radiance (step 2) and then multiply by the solid angle measure of the aperture as seen by a specific point $x_{-1}$ (step 3 in your aperture/pixel reversed scenario has the solid angle of the pixel), and the area of the pixel (step 4 has the area of the aperture). This flip of the aperture and film pixel is not really an issue as noted, but I'll stick to the physical convention (i.e. the reverse of yours). The main issue is that averaging and then multiplying by the solid angle corresponding to the the aperture as seen by the whole pixel, vs by a single point is slightly wrong. Note that $\mathcal{A}_{\mathcal{H}^2(x_{-1}, N_{x_{-1}})}$ is defined per point, and it varies for the different points on the film's pixel, so you cannot precompute one value for the pixel and call it a day, unless the pixel consists only of a single point - then it would be valid.
There is another formulation that is amenable to a slight modification of your idea, however it uses the area formulation instead of the solid angle one. Namely
$$I_{ij} = \int_{t_0}^{t_1}\int_{\mathcal{F}}\int_{\mathcal{A}}w_{ij}(x_{-1}, x_{-1/2})L_i(x_{-1}, -\omega_o))G(x_{-1},x_{-1/2})\, d x_{-1/2}dx_{-1}dt.$$
Here $x_{-1/2}$ is a point on the aperture, and $-\omega_o = \frac{x_{-1/2}-x_{-1}}{\|x_{-1/2}-x_{-1}\|}$, and $G(x_{-1},x_{-1/2}) = \frac{(N_{x_{-1}}\cdot(-\omega_o))(N_{x_{-1/2}}\cdot \omega_o)}{\|x_{-1/2}-x_{-1}\|^2}$. Choose
$$w_{ij}(x_{-1}, x_{-1/2}) = \frac{\boldsymbol{1}_{\mathcal{P}_{ij}}(x_{-1})}{G(x_{-1},x_{-1/2})},$$
where $\mathcal{P}_{ij}\subseteq \mathcal{F}$ is the set of pixel $(i,j)$ on the film, and $\boldsymbol{1}_{\mathcal{X}}$ is the indicator function. Then you get what you want, with the small change that you multiply the average by the areas $|\mathcal{A}|$ and $|\mathcal{P}_{ij}|$ , without involving solid angles of those.
A general remark I have is that if you want to convert between different units then you need integration and differentiation. The multiplication and division is a special case of this, which is valid only in very specific scenarios. My takeaway from your explanation is that this was the main issue with your understanding of the problem.