I am trying to build the mathematical background for an stochastic simulation. Right now I have the mathematical model for:

$$y_t = N(\mu_t, \sigma_t)$$

So that $\mu_t$ and $\sigma_t$ are known. These $y_t$ values are non-independent and are used to evaluate another meassure, let´s say:

$$I_t = (1+y_t) \cdot I_{t-1}$$

Where $I(0)$ is a constant. I am trying to predict the standard deviation of $I_t$ for every $t$.

So far I have tried using the following equations: $$c \cdot N(\mu, \sigma^2) = N(c\mu, (c\sigma)^2)$$ $$c + N(\mu, \sigma^2) = N(c+\mu, \sigma^2)$$ $$N(\mu_1, \sigma_1^2) \cdot N(\mu_2, \sigma_2^2) = N\left(\frac{\sigma_1^2\mu_2+\sigma_2^2\mu_1}{\sigma_1^2+\sigma_2^2}, \frac{1}{\frac{1}{\sigma_1^2}+\frac{1}{\sigma_2^2}}\right)$$

However, after simulating $y_t$ and then converting it to $I_t$ I do not get the same distribution as, with the equations above, I should.

Am I missing something?

EDIT:

I have been able to obtain the mathematical formulation when $\mu_t=0$ $\forall t$ with the help of this link https://mathworld.wolfram.com/NormalProductDistribution.html

To obtain for $\mu_t\neq 0$ I have changed the variables in the integral in the link to: $$P_{XY}(u)=\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}N_x(0,\sigma_x)N_y(0, \sigma_y)\delta\left((x+\mu_x)(y+\mu_y)-u\right)dxdy$$

How does that affect the Meijer G function the obtain?