Advancing a causal IIR filter's impulse response
Let's consider a causal IIR filter with the following transfer function:

$$H(z)=\frac{B(z)}{A(z)}=\frac{\displaystyle\sum_{n=0}^{M}b[n]z^{-n}}{1+\displaystyle\sum_{n=1}^{N}a[n]z^{-n}}\tag{1}$$

Note that $M$ need not be equal to $N$. Let $K$ be the number of samples that must be removed from the impulse response corresponding to the system $(1)$, so the new desired impulse response is

$$\tilde{h}[n]=h[n+K],\qquad n=0,1,\ldots\tag{2}$$

Let $H_K(z)$ denote the FIR transfer function corresponding to the first $K$ samples of $h[n]$:

$$H_K(z)=h[0]+h[1]z^{-1}+\ldots +h[K-1]z^{-(K-1)}\tag{3}$$

The transfer function of the new filter is then given by

$$\tilde{H}(z)=z^K\left[\frac{B(z)}{A(z)}-H_K(z)\right]=\frac{z^K\big[B(z)-H_K(z)A(z)\big]}{A(z)}\tag{4}$$

As expected, the denominator polynomial remains unchanged, we just get a new numerator polynomial

$$\tilde{B}(z)=z^K\big[B(z)-H_K(z)A(z)\big]\tag{5}$$

Let's first consider the case $K=1$. In this case, the FIR transfer function $H_1(z)$ defined in $(3)$ is trivially given by

$$H_1(z)=h[0]=b[0]\tag{6}$$

Consequently, the new numerator polynomial becomes

$$\tilde{B}(z)=\big(b[1]-b[0]a[1]\big)+\big(b[2]-b[0]a[2]\big)z^{-1}+\ldots +\\+\big(b[N]-b[0]a[N]\big)z^{-(N-1)}\tag{7}$$

if $M=N$. If $M&gt;N$, there are additional terms $b[k]z^{-(k-1)}$, $k=N+1,\ldots,M$, and if $M&lt;N$, there are additional terms $-b[0]a[k]z^{-(k-1)}$, $k=M+1,\ldots,N$.

The computation of the new numerator coefficients according to $(7)$ is very simple and efficient. For $K&gt;1$ this procedure can be applied iteratively. In Matlab/Octave this is (almost) a one-liner:


for k = 1:K
    b = [b(2:M+1);zeros(N-M,1)] - b(1)*[a(2:N+1);zeros(M-N,1)];
    M = length(b) - 1;
end


Note that for any $K$, the order of the new numerator polynomial $\tilde{B}(z)$ equals $\max(M-K,N-1)$.


Original Answer (correct, but the presentation and solution given above is probably better; the examples below are still valid and hopefully useful)

A solution can be obtained by recomputing the numerator coefficients as follows:

$$\tilde{b}[n]=(\tilde{h}\star a)[n],\qquad n=0,1,\ldots,L-1\tag{3}$$

where $\star$ denotes convolution, $L=\max(M,N)$, and $\tilde{h}[n]$ is related to the impulse response $h[n]$ of the original IIR filter by

$$\tilde{h}[n]=h[n+1],\qquad n=0,1,\ldots,L-1\tag{4}$$

The method can be summarized as follows:


Compute $L=\max(M,N)$ coefficients of the original impulse response $h[n]$ for indices $n=1,2,\ldots,L$
Compute the new numerator coefficients from the impulse response and from the denominator coefficients using Eq. $(4)$


I chose a random filter to show the result of the suggested procedure in the figure below (blue: original impulse repsonse: red: new impulse response shifted to the right by one sample for comparison):



The left column below shows the original numerator coefficients alongside the modified numerator coefficients, aligned in such a way that it becomes obvious that for higher indices the new coefficients equal the original coefficients. (This is the case because in this example $M&gt;N$; c.f. Eq. $(7)$ of the new answer above, including the text below Eq. $(7)$).



  -0.67010   
   0.65673   0.92477
  -0.61219  -0.74621
  -0.60645  -0.53944
  -0.50316  -0.50316
   0.82855   0.82855
   0.20396   0.20396
   1.15680   1.15680
  -0.29830  -0.29830
   0.10605   0.10605



The denominator coefficients of both filters are


a = [1, 0.4, -0.2, 0.1]



EDIT:

This method can be generalized to remove $K$ initial values of the original impulse response, even if $K&gt;M$:


Compute $L=\max(M,N)$ coefficients of the original impulse response $h[n]$ for indices $n=K,K+1,\ldots,K+L-1$
Compute the new numerator coefficients according to Eq. $(4)$ with $\tilde{h}[n]=h[n+K],\qquad n=0,1,\ldots,L-1$


The resulting numerator coefficients will generally have several trailing zeros (or values sufficiently close to zero) which can be removed.

The following example shows the result of removing $20$ initial values of the impulse response of the filter with coefficients


[b,a] =

   1.000000   1.000000
   2.000000  -0.187178
   3.000000  -0.075960
  -1.000000   0.812250




The numerator coefficients of the new filter are


   0.95236
   0.90498
   0.42252


Of course, these coefficients can no longer be interpreted as (an approximation of) a shifted version of the original numerator coefficients.