Calculate average from a power spectrum
By Welch's method, you can calculate the power spectrum by averaging the magnitude of a bunch of FFT frames (or time frames of the STFT). If your signal is $x(t)$, and its STFT is $X(\omega, \tau)$, where $\omega$ is the frequency bin, $\tau$ is the time frame and $T$ is the total number of frames, the PSD is
$$P_x(\omega) = \frac{1}{T} \sum_{\tau=0}^{T-1} |X(\omega, \tau)|^2 $$
When calculating the STFT, you need to make sure you get constant overlap add. This means, with the Hann window, you need a hopSize of N/2 (where N is your window length). I'd also keep the window length shorter, about 256/512 samples. This is not compulsory in your case, because your signal is stationary. Had it not been stationary you should have chosen a smaller window length. Here is your code, modified.
    y, sr = librosa.load('Noise Recorder.wav',sr=None,duration=2) 
    n_fft = 8192
    win_len = 512
    S = librosa.stft(y=y,n_fft=n_fft,win_length=win_len,hop_length=win_len/2,window='hann')     # Short FFT
    F = librosa.fft_frequencies(sr=sr, n_fft=n_fft) # Frequency components
    Abs = np.abs(S)  # Same as mag, phase = librosa.magphase(S)
    
    PSD = np.mean(Abs**2, axis=1).  #this is the power spectrum
    PSD_dB = 10*np.log10(PSD). #power spectrum in dB
    RMS = np.sqrt(np.mean(PSD**2))) #RMS of the power spectrum
    RMS_dB = 20*np.log10(RMS) # RMS of power spectrum in dB

    ##Theoretical RMS (in dB) for a sine wave of amplitude A dB
    A_lin = 10**(A/20.0)           #convert dB to linear amplitude, A = 114 in your case
    RMS_exp = 20*np.log10(A_lin/n_fft)   #expected RMS in dB
    


`