[scilab-Users] vectorizing for-loops for autoregressive or other recursive processes

[Serge Steer]

Le 21/07/2011 12:18, Ginters Bušs a écrit :
> Dear all,
>
> I am inspired by the recent Allan's comment that I can vectorize 
> for-loops.
>
> Say, I have a for-loop for AR(2) process 
> y(t+1)=f1*y(t)+f2*y(t-1)+epsilon(t+1) with a degenerate disturbance term:
>
> f1=1.2; f2=-0.5;
> n=100000;
> c=1:n;
> c(2)=f1*c(1)+c(2);
> for i=3:n;
> c(i)=c(i)+f1*c(i-1)+f2*c(i-2);
> end
>
> Allan proposed trying
>
> f1=1.2; f2=-0.5;
> n=100000;
> b=1:n;
> b(2)=f1*b(1)+b(2);
> b(3:n) = b(3:n)+f1*b(2:n-1)+f2*b(1:n-2);
>
>
> but the result is different. Then, I thought I can try to use ode; 
> say, for AR(1) process y(t+1) = a*y(t)+ epsilon(t+1) it would look like
>
> deff("yp=a_function(k,y)","yp=a*y+sigma*u(k)")
> y0=0;
> a=.9;
> n=100000;
> u=1:n;
> y=ode("discrete",y0,1,1:n,a_function);
>
> but it turns out that for-loop is a bit faster than this ode code.
>
> Another idea would be using Wold representation of (only stationary) 
> AR process by rewriting the above AR(2) process as
>
> y(t+1) = ((1-f1*L - f2*L^2)^(-1))*epsilon(t+1)
>
> where (1-f1*L - f2*L^2) is a lag polynomial, L defined as Ly(t)=y(t-1)
>
> and trying to get the series representation of (1-f1*L - f2*L^2)^(-1) 
> but I'm stuck here, and this would work only for a stationary process.
>
> Do you have any ideas/experience with rewriting for-loops for AR or 
> other recursive processes more efficiently than the above first code?
>
>
You can use the arsimul Scilab function or create a discrete transfer 
function some thing like

H=syslin('d',1,poly([f2 f1 1],'z','c'))
and the use flts withe ones(1,n) as an imput signal

Serge Steer
INRIA