This can be done, but it would be terribly unwieldy and might be numerically unstable if the predictors are highly correlated.
It's simplest if we assume that each of the $X_i$ and $Y$ values is expressed as its difference from its respective mean. Then we can dispense with the intercept term of the regression equation and the column of 1s more generally included in the design matrix $\mathbf{X}$ (see this page) and write the general solution to the multiple regression in matrix form as:
$$\mathbf{\hat\beta = (X^TX)^{-1}X^T}Y,$$
where $\mathbf{X}$ is the matrix of (mean-centered) predictor values (number of rows $n$ equal to the number of observations, and number of columns $p$ equal to the number of predictors) and $Y$ is the vector of $n$ (mean-centered) outcome observations.
To put this into the context and notation of your question, note that the correlation coefficient $r_{XY}$ is $\mathsf{Cov}(X,Y)/(\sigma_X \sigma_Y),$ where $\mathsf{Cov}(X,Y)$ is the covariance between $X$ and $Y$. For a single predictor, $X$ is a vector with $n$ elements, and with $X$ and $Y$ centered around their mean values
$$X^TY = n \mathsf{Cov}(X,Y).$$
With a single predictor centered around its mean, $\mathbf{X^TX}$ is simply $n \sigma_X^2$, with inverse $1/(n \sigma_X^2)$. By the general matrix formula, $\hat\beta = \mathsf{Cov}(X,Y)/\sigma_X^2$, while your formula gives the same result:
$$ \hat\beta = r_{XY} \frac{\sigma_Y}{\sigma_X} = \frac{\mathsf{Cov}(X,Y)}{\sigma_X \sigma_Y}\frac{\sigma_Y}{\sigma_X}=\frac{\mathsf{Cov}(X,Y)}{\sigma_X^2}.$$
As you increase the number of predictors $p$, $\mathbf{X^T}Y$ is just a simple generalization of the above: it's a vector of the $p$ individual covariances of the $X_i$ with $Y$, multiplied by $n$. If you wish, you can express each covariance in terms of the correlation coefficient and the square roots of individual variances:
$$\mathsf{Cov}(X_i,Y) = r_{X_iY}\sigma_{X_i}\sigma_Y.$$
The inner-product matrix product $\mathbf{X^TX}$ of mean-centered predictor values is a less-simple generalization of the variance of the single predictor; it's $n$ times the variance-covariance matrix of $\mathbf{X}$. That variance-covariance matrix has the individual variances of the $X_i$ along the diagonal, with each off-diagonal $(i,j)$ element equal to the corresponding covariance, $\mathsf{Cov}(X_i,X_j)$. You can express each of these covariances between predictors in terms of the corresponding correlation coefficients and individual (square roots of) variances, similarly to the covariance between an individual predictor and the outcome as above.
You must, however, take the inverse of the predictor variance-covariance matrix to calculate the vector of regression coefficients $\mathbf{\beta}.$ That's where expressing things in terms of correlation coefficients and individual (square roots of) variances gets complicated.
You could get a closed-form solution for $\mathbf{X^TX}$ for any number of predictors with elements expressed in that way. Brute force with the formula for a matrix inverse would do the trick, or you could start with $\mathbf{(X^TX)^{-1}}=1/(n \sigma_x^2)$ for the single-predictor case and add one predictor at a time following the formula on Page 19 of the Matrix Cookbook for updating $\mathbf{(X^TX)^{-1}}$ when you add one predictor vector to it. This gets messy very quickly, however.
Even if you got a closed-form solution in terms of correlation coefficients and individual variances, it would be unwise just to plug into that formula particularly if there were strong correlations among the predictors. In that case you might get numerically unstable results. Standard statistical software is designed to solve this matrix equation in a way that minimizes those problems.
Heuristically, think of the multiple-predictor formula
$$\mathbf{\hat\beta = (X^TX)^{-1}X^T}Y$$
as a generalization of the single-predictor formula
$$ \hat\beta =\frac{\mathsf{Cov}(X,Y)}{\sigma_x^2}.$$
In both cases, you are multiplying the covariance(s) between the outcome and the predictor(s) by the inverse of the (co)variance(s) of the predictors. That's perhaps simpler than thinking in terms of correlation coefficients.* To do the calculation in practice though, let well vetted statistical software to the work for you.
In response to update and comment: As you are working with time-series data this approach is probably inappropriate for your application. That's particularly true if you are estimating the variances and the correlation coefficients from different data windows.
*If you want to think in terms of correlation coefficients, start with outcomes and predictors each not only centered around their means but also scaled to unit standard deviation/variance. Then the diagonal entries of $\mathbf{X^TX}$ are all 1, the off-diagonal entries are the pairwise correlation coefficients between the predictors, and $\mathbf{X^T}Y$ is the vector of correlation coefficients between each of the predictors and $Y.$ After applying the formula for $\mathbf{\hat\beta}$ re-scale the regression coefficients to take the initial scaling into account; that’s where the factors like $\sigma_Y/\sigma_{X_i}$ come from in your formula for the two-predictor case. You could use a symbolic mathematics manipulation program like Mathematica to do the algebra. I gave that a try, recovering your formula for the 2-predictor case. The result wasn't terribly complicated for 3 predictors, but with 4 or 5 predictors you end up with very many terms both in the numerator and denominator.
