You try to make a stochastic model corresponding to the deterministic model above. The Master equation of that model is
$$\begin{eqnarray}\frac{d}{dt}\mathbb{P}(I,P) & = & \mathbb{T}(I,P|I-1,P)\cdot \mathbb{P}(I-1,P) - \mathbb{T}(I+1,P|I,P)\cdot \mathbb{P}(I,P) \\
& & + \mathbb{T}(I,P|I+1,P)\cdot \mathbb{P}(I+1,P) - \mathbb{T}(I-1,P|I,P)\cdot \mathbb{P}(I,P) \\
& & + \mathbb{T}(I,P|I,P-1)\cdot \mathbb{P}(I,P-1) - \mathbb{T}(I,P+1|I,P)\cdot \mathbb{P}(I,P) \\
& & + \mathbb{T}(I,P|I,P+1)\cdot \mathbb{P}(I,P+1) - \mathbb{T}(I,P-1|I,P)\cdot \mathbb{P}(I,P) \; .
\end{eqnarray}$$
Where the transition probabilities are
$$\begin{eqnarray}
\mathbb{T}(I+1,P|I,P) & = & \alpha_I+P\lambda_I\\
\mathbb{T}(I-1,P|I,P) & = & I\gamma_I-PI\delta_I\\
\mathbb{T}(I,P+1|I,P) & = & \alpha_P+P\theta_P\\
\mathbb{T}(I,P-1|I,P) & = & PI\delta_P
\end{eqnarray}$$
It takes a bit of work, but computing the expectation values $\mathbb{E}(I)=\sum_{I,P} I\cdot \mathbb{P}(I,P)$ and $\mathbb{E}(P)=\sum_{I,P} P \cdot\mathbb{P}(I,P)$, can be done by multiplying the master equation by either $I$ or $P$ and then summing over all values that $I$ and $P$ can get. This will give you the following equations:
$$\begin{cases}
\frac{\mathrm{d}\mathbb{E}(I)}{\mathrm{d}t}=& \alpha_I+\lambda_I\mathbb{E}(P)- \gamma_I\mathbb{E}(I) -\delta_I\mathbb{E}(PI) \\
\frac{\mathrm{d}\mathbb{E}(P)}{\mathrm{d}t}=&\alpha_P+\theta_P\mathbb{E}(P)-\delta_P\mathbb{E}(PI)
\end{cases}$$
Note however that you get the expectation of the product $PI$ and not the product of expectations.
EDIT:
Here's my code in R for a run of the stochastic model
ai<-0; # parameters
gi<-0.1;
di<-0.05;
li<-0.05;
tp<-0.15;
dp<-0.01;
ap<-0.5;
T<-10;
plusI<-c(1,-1,0,0);
plusP<-c(0,0,1,-1);
Time<-c(0); # initial conditions
I<-c(20);
P<-c(5);
for (k in 1:800) {
freq<-ai+li*P[k]+gi*I[k]+di*P[k]*I[k]+ap+tp*P[k]+dp*P[k]*I[k];
t<-rexp(1,freq);
Time<-c(Time,Time[k]+t);
sel<-sample(c(1,2,3,4),1,prob=c(ai+li*P[k],gi*I[k]+di*P[k]*I[k],ap*(Time[k]
To explain the most essential part, I draw a time from an exponential distribution with parameter the sum of all the rates of the processes. Then, I decide which process takes place by sampling with a probability vector equal to the 4 rates of the processes. This way, I manage to simulate the stochastic process. Here's a run of it:

I haven't tried to explore the state space in detail yet, but all choices of initial conditions I have tried lead to an exponential explosion of the rabbit population and a collapse of the fox population.