Convolution to do FIR with initial and final conditions!
Your benchmark is actually quite surprising! I would've expected filter to choose an appropriate implementation of convolution to filter the input data, but thinking about this, since filter needs to pass the state in a format that isn't dependent on the convolution implementation actually used, the whole process of generating the state might take longer than the actual filtering; note that this really is just speculation. (Assumption: filter with long input data uses FFT FIR implementation internally, but the state that filter returns are taps of a direct-form FIR implementation – you basically need to reverse quite a bit of the filtering you've done to get that)


  I already tried to overlap the chunks of data and mimic the "state" for the convolution but haven't been successful yet.


What you probably would like to have is a well-implemented FFT FIR! Since b doesn't seem to change, this would be a very efficient way of implementing the filtering operation.

The basic idea would be described as overlap-add filter (wikipedia), and the rough idea is this:


Pad your tap vector $b$ with zeros to length $N$. Get the $N$-point discrete Fourier transform of that using the FFT. Call the result $H$.
Take segments of length $L$ of the input signal, and pad them to $N$, transform as above.
point-wise multiply, since multiplication in discrete freq domain is equivalent to circular convolution in time domain.
inverse DFT (via IFFT)
By zero-padding enough, we now that a $L$ long part of the output signal can be generated by adding up the $N-L$ last samples of the last iteration to the start of the samples of this iteration.
save the $N-L$ last output samples for next iteration


 

personal addition: If you already know what filter you want, and you're actually after speed and don't care as much that you might need to leave Matlab to do the filtering (I personally cherish every moment that I don't have to use Matlab), look at this explanatory blog post that I've written, which actually does a lot of of band pass filters in parallel – you'd only need one of these, and you could choose a FFT FIR block to do the work; I've been doing 160 megasamples per second with these, on a "normal" PC. This is, as I'll frankly admit, based on GNU Radio and the VOLK acceleration library, because that's what I prefer for streaming samples DSP (which you seem to be doing in Matlab).