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

[Ginters Bušs]

On Thu, Jul 21, 2011 at 3:21 PM, Serge Steer <Serge.Steer at inria.fr> wrote:

> 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
>
>
>
> ok, arsimul and armac look quite bulky for vector size 1e5, and the help is
not that helpful without showing the model; it is quite difficult to
understand what is what without seeing the model;

A different thing is about syslin - it appears I like this stuff, but I can
not get the wanted result, may be because I do not understand what I am
doing.

Particularly, if I write (a stationary) AR(2) model as x_t =((1-f1*L -
f2*L^2)^(-1))*epsilon_t, then the polynomial is not
H=syslin('d',1,poly([f2 f1 1],'z','c'))
but is
H=syslin('d',1,poly([1 -f1 -f2],'z','c')).

Also, in the example above, I used an input signal as c=1:n, not
c=ones(1,n), so shouldn't I use
y=flts(1:n,H);
not
y=flts(ones(1,n),H);?

However, trying various combinations of the former and the latter, I cannot
replicate this:

b=1:n;
b(2)=f1*b(1)+b(2);
for i=3:n;
    b(i)=b(i)+f1*b(i-1)+f2*b(i-2);
end

may be this is because of initial undeclared conditions b_0=0; b_{-1}=0 in
the flts command; therefore, I tried to compare it with

b=1:n;
b(1)=0;
b(2)=0;
for i=3:n;
    b(i)=b(i)+f1*b(i-1)+f2*b(i-2);
end

Still different. No wonder this is the case because have little clue what I
am doing.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://lists.scilab.org/pipermail/users/attachments/20110721/24d079e7/attachment.htm>