Output power of a signal passing through an IIR filter suppose I have the following discrete time signal
$x\left(n\right)=0.7\sin\left(2\pi \times 623\:nT\right)$
where $N=4000$ and $T=\frac{1}{f_s}$, where $f_s=40$ kHz
I pass this signal through a simple IIR filter $ H\left(z\right)\:=\:\frac{0.2}{1-0.3z^{-1}}$
I calculate the power of the ouput signal $P_y\:=\:P_x\:\left|H\left(z\right)\right|^2_2$
basically am using the filter 2-norm in this case.
for that my calculations would be as follows:
$P_y\:=\left(0.245\right)\left(\frac{0.2^2}{1-0.3^2}\right)\:=\:0.0108$
However when I develop this in Matlab, I get the output power as 0.0199
Here's my code in Matlab
fs = 40000;
f1 = 623;
N = 4000;
n = 1:N;
b = 0.2;
a = [1 -0.3];
x = 0.7*sin(2*pi*f1*n/fs);
P_x = var(x);
y = filter(b,a,x);
P_y = var(y);

I don't understand what I was doing wrong, would someone please clarify for me? I would appreciate it. Thank you in advance.
The analog frequency $\Omega=2\pi\times 623$ and the digital frequency is $\omega = \Omega T = \Omega/f_s = 0.03115\pi$
The frequency response at this frequency is
$$H(z)|_{z=e^{j\omega}} = H(e^{j0.03115\pi}) = \frac{0.2}{1-0.3 e^{-j0.03115\pi}}$$
and its square of magnitude is
$$
|H(e^{j0.03115\pi})|^2 = \Big|\frac{0.2}{1-0.3 e^{-j0.03115\pi}}\Big|^2 = 0.0812
$$
Now you'll see $P_y = P_x |H(e^{j0.03115\pi})|^2 = 0.0199$