new file mode 100644
--- /dev/null
+++ b/scripts/signal/arch_rnd.m
@@ -0,0 +1,87 @@
+## Copyright (C) 1995, 1996, 1997  Kurt Hornik
+## 
+## This program is free software; you can redistribute it and/or modify
+## it under the terms of the GNU General Public License as published by
+## the Free Software Foundation; either version 2, or (at your option)
+## any later version.
+## 
+## This program is distributed in the hope that it will be useful, but
+## WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+## General Public License for more details. 
+## 
+## You should have received a copy of the GNU General Public License
+## along with this file.  If not, write to the Free Software Foundation,
+## 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
+
+## usage:  y = arch_rnd (a, b, T)  
+##
+## Simulates an ARCH sequence y of length T with AR coefficients b and
+## CH coefficients a.
+## I.e., y follows the model
+##     y(t) = b(1) + b(2) * y(t-1) + ... + b(lb) * y(t-lb+1) + e(t),
+## where e(t), given y up to time t-1, is N(0, h(t)), with
+##     h(t) = a(1) + a(2) * e(t-1)^2 + ... + a(la) * e(t-la+1)^2.
+
+## Author:  KH <Kurt.Hornik@ci.tuwien.ac.at>
+## Description:  Simulate an ARCH process
+
+function y = arch_rnd (a, b, T)
+  
+  if (nargin != 3)
+    usage ("arch_rnd (a, b, T)");
+  endif
+  
+  if !( (min (size (a)) == 1) && (min (size (b)) == 1) )
+    error ("arch_rnd:  a and b must both be scalars or vectors");
+  endif
+  if !( is_scalar (T) && (T > 0) && (rem (T, 1) == 0) )
+    error ("arch_rnd:  T must be a positive integer");
+  endif
+  
+  if !(a(1) > 0)
+    error ("arch_rnd:  a(1) must be positive");
+  endif
+  ## perhaps add a test for the roots of a(z) here ...
+  
+  la = length (a);
+  a  = reshape (a, 1, la);
+  if (la == 1)
+    a  = [a, 0];
+    la = la + 1;
+  endif
+  lb = length (b);
+  b  = reshape (b, 1, lb);
+  if (lb == 1)
+    b  = [b, 0];
+    lb = lb + 1;
+  endif
+  M  = max([la lb]);
+  
+  e  = zeros (T, 1);
+  h  = zeros (T, 1);
+  y  = zeros (T, 1);
+  
+  h(1) = a(1);
+  e(1) = sqrt (h(1)) * randn;
+  y(1) = b(1) + e(1);
+  
+  for t= 2 : M;
+    ta   = min ([t la]);
+    h(t) = a(1) + a(2:ta) * e(t-1:t-ta+1).^2;
+    e(t) = sqrt (h(t)) * randn;
+    tb   = min ([t lb]);
+    y(t) = b(1) + b(2:tb) * y(t-1:t-tb+1) + e(t);
+  endfor
+  if (T > M)
+    for t = M+1 : T;
+      h(t) = a(1) + a(2:la) * e(t-1:t-la+1).^2;
+      e(t) = sqrt (h(t)) * randn;
+      y(t) = b(1) + b(2:lb) * y(t-1:t-tb+1) + e(t);
+    endfor
+  endif
+  
+  y = y(1:T);
+  
+endfunction
+    
\ No newline at end of file