function [rsim,xsim] = DGPsimuZero(cvec,gvec,T,Sigma,Nsimu)

rho1 = 1-cvec(1)/T;
rho2 = 1-cvec(2)/T;
theta = 1-gvec/T;

[arcoef,macoef] = root2coef(rho1,rho2,theta);

arcoef = -arcoef;

%% simulations %%
errsim = mvnrnd(zeros(2,1),Sigma,T*Nsimu);

rsim = vec2mat(errsim(:,1)',T);
rsim = rsim';

esim = vec2mat(errsim(:,2)',T);
esim = esim';

if rho1 == theta
    xsim = filter(1,[1 -rho2],esim,[],1);
elseif rho2 == theta 
    xsim = filter(1,[1 -rho1],esim,[],1);
else
    xsim = filter([1 macoef],[1 arcoef],esim,[],1);
end


end










