IIR implementation - RuntimeWarning: overflow encountered in double_scalars
First, I assume that it is a that is equal to array([ 1.        , -3.9671626 ,  5.90202586, -3.90255878,  0.96769554]).
Second, that array when you go through all the z-domain stuff, you find that you are representing the difference equation
$$\sum_{k=0}^N a_k y_{n-k} = \sum_{k=0}^N b_k x_{n-k}$$
If you solve that for $y_n$ you get
$$y_n = -\sum_{k=1}^N a_k y_{n-k} + \sum_{k=0}^N b_k x_{n-k}$$
Note the minus sign in front of the summation of the prior terms of $y$ -- your code is missing that.  So you've turned your denominator polynomial from $a(z) = 1 - 3.967z^{-1} + 5.902z^{-2} \cdots$ into $a(z) = 1 + 3.967z^{-1} - 5.902z^{-2} \cdots$.  That polynomial is almost certainly unstable, and you're overflowing even what double-precision floating point can represent.
Note 1:
In general, it's better to break filters down into the smallest sections you can and cascade them.  This is because (a) it is the filter poles and zeros that are important to filter behavior, but polynomial roots get absurdly sensitive to polynomial coefficient variation as the order goes up, and (b) the filter's susceptibility to rounding in the data -- even with perfect coefficients -- goes up pretty much the same as a polynomial's sensitivity to coefficient variation.
So -- small polynomials are good, which means low-order filters are good.
Your filter design function should be able to cough up a form where it gives you poles &amp; zeros instead of polynomials; you can then break those into 2nd-order functions.
Note 2:
You're using Scipy, which -- if I'm not mistaken -- has a built-in filtering function, and it should be compatible with the poles &amp; zeros form of the filter design output.  If you're rolling your own filter code to learn it -- good for you, keep it up, we need people to know what's going on under the hood.  If you just need to filter some @#$% data and get on with life, use the built-in Scipy function.