<br><br><div class="gmail_quote">On Thu, Jul 21, 2011 at 3:21 PM, Serge Steer <span dir="ltr"><<a href="mailto:Serge.Steer@inria.fr">Serge.Steer@inria.fr</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex;">
Le 21/07/2011 12:18, Ginters Bušs a écrit :<br>
<div><div></div><div class="h5">> Dear all,<br>
><br>
> I am inspired by the recent Allan's comment that I can vectorize<br>
> for-loops.<br>
><br>
> Say, I have a for-loop for AR(2) process<br>
> y(t+1)=f1*y(t)+f2*y(t-1)+epsilon(t+1) with a degenerate disturbance term:<br>
><br>
> f1=1.2; f2=-0.5;<br>
> n=100000;<br>
> c=1:n;<br>
> c(2)=f1*c(1)+c(2);<br>
> for i=3:n;<br>
> c(i)=c(i)+f1*c(i-1)+f2*c(i-2);<br>
> end<br>
><br>
> Allan proposed trying<br>
><br>
> f1=1.2; f2=-0.5;<br>
> n=100000;<br>
> b=1:n;<br>
> b(2)=f1*b(1)+b(2);<br>
> b(3:n) = b(3:n)+f1*b(2:n-1)+f2*b(1:n-2);<br>
><br>
><br>
> but the result is different. Then, I thought I can try to use ode;<br>
> say, for AR(1) process y(t+1) = a*y(t)+ epsilon(t+1) it would look like<br>
><br>
> deff("yp=a_function(k,y)","yp=a*y+sigma*u(k)")<br>
> y0=0;<br>
> a=.9;<br>
> n=100000;<br>
> u=1:n;<br>
> y=ode("discrete",y0,1,1:n,a_function);<br>
><br>
> but it turns out that for-loop is a bit faster than this ode code.<br>
><br>
> Another idea would be using Wold representation of (only stationary)<br>
> AR process by rewriting the above AR(2) process as<br>
><br>
> y(t+1) = ((1-f1*L - f2*L^2)^(-1))*epsilon(t+1)<br>
><br>
> where (1-f1*L - f2*L^2) is a lag polynomial, L defined as Ly(t)=y(t-1)<br>
><br>
> and trying to get the series representation of (1-f1*L - f2*L^2)^(-1)<br>
> but I'm stuck here, and this would work only for a stationary process.<br>
><br>
> Do you have any ideas/experience with rewriting for-loops for AR or<br>
> other recursive processes more efficiently than the above first code?<br>
><br>
><br>
</div></div>You can use the arsimul Scilab function or create a discrete transfer<br>
function some thing like<br>
<br>
H=syslin('d',1,poly([f2 f1 1],'z','c'))<br>
and the use flts withe ones(1,n) as an imput signal<br>
<font color="#888888"><br>
Serge Steer<br>
INRIA<br>
<br>
<br>
<br>
</font></blockquote></div>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;<br>
<br>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. <br><br>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 <br>
H=syslin('d',1,poly([f2 f1 1],'z','c'))<br>but is <br>H=syslin('d',1,poly([1 -f1 -f2],'z','c')).<br><br>Also, in the example above, I used an input signal as c=1:n, not c=ones(1,n), so shouldn't I use <br>
y=flts(1:n,H);<br>not<br>y=flts(ones(1,n),H);? <br><br>However, trying various combinations of the former and the latter, I cannot replicate this:<br><br>b=1:n;<br>b(2)=f1*b(1)+b(2);<br>for i=3:n;<br>    b(i)=b(i)+f1*b(i-1)+f2*b(i-2);<br>
end<br><br>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 <br><br>b=1:n;<br>b(1)=0;<br>b(2)=0;<br>for i=3:n;<br>    b(i)=b(i)+f1*b(i-1)+f2*b(i-2);<br>
end <br><br>Still different. No wonder this is the case because have little clue what I am doing.<br>