Humour me please. Let us suppose we have an ARMAX(1,0,0) model with 2 independent variables and a constant.

ie  (1-ΦB)Yt = α +µ1*X1t +µ2*X2t  where B is the lag or backshift operator

Assuming Φ is positive this is equivalent to 

Yt =α + Φ*Yt-1 +µ1*X1t +µ2*X2t 

  which is relatively straightforward to check against Gretl output  using the estimated coefficients and relevant rhs variables .However, the results of doing so do not correspond to the predicted or fitted values  from Gretl.

I think section 31.2 of the user guide could be relevant here, which talks explicitly about ARMAX models and their subtleties.

cheers

sven