How do I run MP2 and/or CCSD-level hyperpolarizability calculations with GAMESS?
As you guessed, to calculate $\beta$ within GAMESS for MP2 and CCSD, you would need to perform first- and second-order finite difference, respectively, as there is no coupled cluster gradient needed, so no true CCSD dipole moment.

However, there is a more fundamental problem: it is impossible to calculate time-dependent response using a finite field approach. Take the general form of 1-dimensional 2nd-order central difference,
$$
f''(x) \approx \lim_{h \to 0} \frac{f(x+h)-2f(x)+f(x-h)}{h^2},
$$

which can directly translate to a diagonal element of the polarizability tensor formed from energies at finite electric fields:

$$
\alpha_{zz}(x) \approx \lim_{h \to 0} \frac{E(x+h)-2E(x)+E(x-h)}{h^2}
$$

where $h$ is a finite electric dipole field of some strength applied along the $z$-direction. What is $x$? Is it as simple as converting the desired frequency $\omega$ to units of electric field strength? No; that still corresponds to a static field rather than a time-varying field. How are the different perturbation directions handled? $x$ must be set to 0.

From a molecular properties review by Jurgen Gauss, page 5:


  (...) disadvantages of the numerical differentiation scheme are
  
  a) that there is no straightforward extension to the computation of frequency-dependent properties (...)


Page 37:


  While analytic derivative theory is sufficient for the theoretical treatment of time-independent (static) properties, the underlying theory needs to be extended for the calculation of time-dependent (dynamical) properties. In particular, the fact that there is -- unlike for the static case -- in the time-dependent case no well-defined energy explains why the simple derivative theory discussed so far is not applicable.


Page 39 (with some manipulation by me):


  It can be shown that the linear, quadratic, etc. response functions $\left&lt;\left&lt;A;B\right&gt;\right&gt;_{w_{b}}$, $\left&lt;\left&lt;A;B,C\right&gt;\right&gt;_{w_{b},w_{c}}$, etc. and thus the frequency-dependent properties of interest can be determined as derivatives of the so-called time-averaged quasi energy
  
  $$
Q(t) = \left&lt;\tilde{\psi}\left|\left(H - i\frac{\partial}{\partial t}\right)\right|\tilde{\psi}\right&gt;
$$
  
  with the phase-isolated wavefunction
  
  $$
\tilde{\psi}(t) = e^{iF(t)} \psi(t)
$$


where $F$ is a time-varying external field, and $\psi(t)$ is the usual wavefunction where MO coefficients and thus the density are explicitly time-dependent.

The only numerical (rather than analytic) method for calculating dynamic properties I am aware of is direct time-propagation of the wavefunction. If you want to calculate frequency-dependent CCSD hyperpolarizabilities, your best (free) bet is DALTON. Note that their coupled cluster module is closed-shell only.