Here are a few things to consider.

(1) Actually your formula had a sign error (which i found by comparing HermiteH[4,x] behavior with HermiteH[n,x] and then substituting n->4).

http://functions.wolfram.com/Polynomials/HermiteH/17/01/01/

So the corrected version of my replacement rule code would be

rul = (x_^p_.* HermiteH[n_, x_] /; IntegerQ[p] && p >= 1) :> 
   x^(p - 1)*HermiteH[n + 1, x]/2 + n *x^(p - 1)*HermiteH[n - 1, x];

(2) I do not know how to set up and solve an appropriate partial difference equation to obtain a general form. I suspect there is no closed form since it would involve a sum that likely has no closed form.

(3) An alternative to rule replacement is as follows (works for specific p but general n).

Clear[f]
f[0, n_] := HermiteH[n, x]
f[p_, n_] := f[p, n] = Expand[f[p - 1, n + 1]/2 + n*f[p - 1, n - 1]]

Example:

f[3, n]

(* Out[48]= 2 n HermiteH[-3 + n, x] - 3 n^2 HermiteH[-3 + n, x] + 
 n^3 HermiteH[-3 + n, x] + 3/2 n^2 HermiteH[-1 + n, x] + 
 3/4 HermiteH[1 + n, x] + 3/4 n HermiteH[1 + n, x] + 
 1/8 HermiteH[3 + n, x] *)

Are you looking for a symbolic Sum result for unspecified parameter values of p, or for code that does the decomposition into an explicit Plus expression when given specific integer values for p?

If the latter, could use a replacement rule, applied repeatedly.

rul = (x_^p_.* HermiteH[n_, x_] /; IntegerQ[p] && p >= 1) :> 
   x^(p - 1)*HermiteH[n + 1, x]/2 - n *x^(p - 1)*HermiteH[n - 1, x];

Example:

In[187]:= Expand[x^3*HermiteH[n, x] //. rul]

(* Out[187]= -2 n HermiteH[-3 + n, x] + 3 n^2 HermiteH[-3 + n, x] - 
 n^3 HermiteH[-3 + n, x] + 3/2 n^2 HermiteH[-1 + n, x] - 
 3/4 HermiteH[1 + n, x] - 3/4 n HermiteH[1 + n, x] + 
 1/8 HermiteH[3 + n, x] *)

I'm not interested in calculating the quantum mechanical expectation value, which is the moment of the square of the wave function. I'm using the QM Harmonic oscillator wave functions as a basis set in which to expand a function I want to compute the moments of, so I want the moments of the HO wave functions, not its square.

Moments involve the square of the wave function. For example, the m^th moment is given by Integrate[x^m psi_n(x)^2, {x, -inf, inf}].

Many results involving Hermite polynomials can be found by making use of the generating function exp(2xt-t^2) = sum ...

In my Mathematica Integrate command I included the restriction that n and p both be positive integers.

Since HermiteH[n,x] is a polynomial of degree n, the coefficient of Exp[-x^2/2] in the integrand is a polynomial of degree n+p. For large x, this behaves as x^(n+p) Exp[-x^2/2], which clearly goes to 0 faster than 1/x, so the integral converges..