We propose a mechanism which can improve the numerical robustness of a subspace based system identification method, the PI scheme [14, Part III], when the unknown system has poles situated close to z = 1, a condition that often arises in applications where the sampling rate is too high. The PI method is capable of solving a deterministic MIMO identification problem in which the output can be corrupted by a very general perturbation including arbitrarily colored noise, transients due to nonzero initial conditions, and a deterministic zero bias, By performing a bilinear transformation on the shift operator we are able to move the poles away from the point z = 1 and a more robust identification results. The implementation of this transformation gives rise to a series of anticausal filters applied to the input/output data. Estimation accuracy is further improved by taking the unknown end conditions of the anticausal filters into account, particularly when only short data records are available, A numerical simulation highlights the improvements realized by our new algorithms.