1 files changed, 95 insertions, 22 deletions
diff --git a/macros/arch_rnd.sci b/macros/arch_rnd.sci
index 50facc4..9e3d88c 100644
--- a/macros/arch_rnd.sci
+++ b/macros/arch_rnd.sci
@@ -1,11 +1,13 @@
-function res = arch_rnd (a, b, t)
 //Simulate an ARCH sequence of length t with AR coefficients b and CH coefficients a.
+
 //Calling Sequence
 //arch_rnd (a, b, t)
-//Parameters 
+
+//Parameters
 //a: CH coefficients
 //b: AR coefficients
 //t: Length of ARCH sequence
+
 //Description
 //This is an Octave function.
 //It Simulates an ARCH sequence of length t with AR coefficients b and CH coefficients a.
@@ -15,29 +17,100 @@ function res = arch_rnd (a, b, 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
+
 //Examples
 //a = [1 2 3 4 5];
 //b = [7 8 9 10];
+//t = 5 ;
 //arch_rnd (a, b, t)
-//ans =
+//Output
+// ans  =
+//
+//    7.2113249
+//    65.479684
+//    654.00814
+//    7194.6572
+//    78364.905
+
+
+
+//function res = arch_rnd (a, b, t)
+//funcprot(0);
+//lhs = argn(1)
+//rhs = argn(2)
+//if (rhs < 3 | rhs > 3)
+//error("Wrong number of input arguments.")
+//end
+//
+//select(rhs)
 //
-//   6.1037e+00
-//   5.7294e+01
-//   5.7390e+02
-//   6.3063e+03
-//   6.8695e+04
-
-funcprot(0);
-lhs = argn(1)
-rhs = argn(2)
-if (rhs < 3 | rhs > 3)
-error("Wrong number of input arguments.")
-end
-
-select(rhs)
-	
-	case 3 then
-	res = callOctave("arch_rnd",a, b, t)
-
-	end
+//	case 3 then
+//	res = callOctave("arch_rnd",a, b, t)
+//
+//	end
+//endfunction
+
+
+function y = arch_rnd (a, b, t)
+
+  funcprot(0);
+  [nargout, nargin] = argn() ;
+
+  if (nargin ~= 3)
+    error("arch_rnd: invalid input arguments")
+  end
+
+  if (~ ((min (size (a)) == 1) & (min (size (b)) == 1)))
+    error ("arch_rnd: A and B must both be scalars or vectors");
+  end
+  if (~ (isscalar (t) & (t > 0) & (modulo(t, 1) == 0)))
+    error ("arch_rnd: T must be a positive integer");
+  end
+
+  if (~ (a(1) > 0))
+    error ("arch_rnd: A(1) must be positive");
+  end
+  // perhaps add a test for the roots of a(z) here ...
+
+  la = length (a);
+  a  = matrix(a, 1, la);
+  if (la == 1)
+    a  = [a, 0];
+    la = la + 1;
+  end
+
+  lb = length (b);
+  b  = matrix(b, 1, lb);
+  if (lb == 1)
+    b  = [b, 0];
+    lb = lb + 1;
+  end
+  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)) * rand();
+  y(1) = b(1) + e(1);
+
+  for t = 2:m
+    ta   = min ([t, la]);
+    h(t) = a(1) + a(2:ta) * e(t-ta+1:t-1).^2;
+    e(t) = sqrt (h(t)) * rand() ;
+    tb   = min ([t, lb]);
+    y(t) = b(1) + b(2:tb) * y(t-tb+1:t-1) + e(t);
+  end
+
+  if (t > m)
+    for t = m+1:t
+      h(t) = a(1) + a(2:la) * e(t-la+1:t-1).^2;
+      e(t) = sqrt (h(t)) * rand() ;
+      y(t) = b(1) + b(2:lb) * y(t-tb+1:t-1) + e(t);
+    end
+  end
+
+  y = y(1:t);
+
 endfunction