Is there a similar formula like Ramanunjan's Eisenstein series identity for $\sum_{k=1}^{n-1}k^2 \sigma(k)\sigma(n-k)$?

All these identities can indeed be proved essentially trivially using modular forms and quasi-modular forms (those involving $E_2$), and the fact that the dimension of such spaces is $1$ for weight 4,6,8,10,14, and $2$ for weight 12, in which case the identities involve also the Ramanujan $\tau$ function. Explicitly, sums $\sum_{1\le k\le n-1}k^a\sigma_b(k)\sigma_c(n-k)$ with $b$ and $c$ odd positive integers ($\sigma_b(k)=\sum_{d\mid k}d^b$) have weight $w=b+c+2+2a$, so if $w=4,6,8,10,14$ you will obtain identities involving only $\sigma$, and if $w=12$ also $\tau(n)$.

Numerical experiments suggest that $$A_2(n) := \sum_{k=1}^{n-1} k^2\sigma(k)\sigma(n-k) = \frac{n^2}{8}\sigma_3(n) - \frac{4n^3-n^2}{24}\sigma(n).$$ PS. In fact, it directly follows from the quoted Touchard and Ramanujan identities.

A couple of similar identities: $$A_1(n):=\sum_{k=1}^{n-1} k\sigma(k)\sigma(n-k) = \frac{5n}{24}\sigma_3(n) - \frac{6n^2-n}{24}\sigma(n).$$ $$A_3(n):=\sum_{k=1}^{n-1} k^3\sigma(k)\sigma(n-k) = \frac{n^3}{12}\sigma_3(n) - \frac{3n^4-n^3}{24}\sigma(n).$$

---

ADDED. A recurrent formula for $A_d(n)$ with an odd $d$ can be obtained from the observation: \begin{split} A_d(n) & := \sum_{k=1}^{n-1} k^d\sigma(k)\sigma(n-k) \\ &= \sum_{k=1}^{n-1} (n-k)^d\sigma(k)\sigma(n-k) \\ &= \sum_{i=0}^d \binom{d}{i} n^{d-i} (-1)^i A_i(n). \end{split} implying that \begin{split} A_d(n) &= \frac{1}{2} \sum_{i=0}^{d-1} \binom{d}{i} n^{d-i} (-1)^i A_i(n) \\ &=\frac{1}{d+1} \sum_{i=0}^{d-1} \binom{d+1}{i} n^{d-i} (-1)^i A_i(n). \end{split} However, to use this formula one would need to compute $A_t(n)$ for even $t<d$ by other means.

It also follows that the generating function: $$\mathcal{A}_n(x) := \sum_{d=0}^{\infty} \frac{A_d(n)}{n^d}x^d$$ satisfied the functional equation: $$\mathcal{A}_n(x) = \frac{1}{1-x}\mathcal{A}_n(\frac{x}{x-1}).$$