Merge pull request #17 from avinashlalotra/masterHEADmaster

Abinash's Work

1 files changed, 85 insertions, 35 deletions

diff --git a/macros/arch_fit.sci b/macros/arch_fit.sci
index 3e431a6..2fa462c 100644
--- a/macros/arch_fit.sci
+++ b/macros/arch_fit.sci
@@ -1,35 +1,85 @@
-function [A, B] = arch_fit(Y, varargin)
-//This functions fits an ARCH regression model to the time series Y using the scoring algorithm in Engle's original ARCH paper.
-//Calling Sequence
-//[A, B] = arch_fit (Y, X, P, ITER, GAMMA, A0, B0)
-//Parameters
-//Description
-//Fit an ARCH regression model to the time series Y using the scoring algorithm in Engle's original ARCH paper.
-//
-//The model is
-//
-//          y(t) = b(1) * x(t,1) + ... + b(k) * x(t,k) + e(t),
-//          h(t) = a(1) + a(2) * e(t-1)^2 + ... + a(p+1) * e(t-p)^2
-//
-//in which e(t) is N(0, h(t)), given a time-series vector Y up to time t-1 and a matrix of (ordinary) regressors X up to t. The order of the regression of the residual variance is specified by P.
-//
-//If invoked as 'arch_fit (Y, K, P)' with a positive integer K, fit an ARCH(K, P) process, i.e., do the above with the t-th row of X given by
-//
-//          [1, y(t-1), ..., y(t-k)]
-//
-//Optionally, one can specify the number of iterations ITER, the updating factor GAMMA, and initial values a0 and b0 for the scoring algorithm.
-funcprot(0);
-rhs = argn(2);
-lhs=argn(1);
-if(rhs<7 | rhs>7)
-error("Wrong number of input arguments.");
-end
-if (lhs<2 | lhs>2)
-    error("Wrong number of output arguments.");
-end
-
-	select(rhs)
-	case 7 then
-	[A, B] = callOctave("arch_fit",Y, varargin(1), varargin(2), varargin(3), varargin(4), varargin(5), varargin(6));
-	end
-endfunction
+/*
+Dependencies : ols, autoreg_matrix
+Calling Sequence
+      [a, b] = arch_fit (y, x, p)
+      [a, b] = arch_fit (y, x, p, iter, gamma, a0, b0)
+Parameters
+      y(vector) : A time-series data vector up to time t-1 .
+      x (Matrix): A matrix of (ordinary) regressors x up to t.
+      p (scalar): The order of the regression of the residual variance.
+      iter (scaler) : Number of iterations
+      gamma (real number) : updating factor
+      a0 b0 (real numbers) : Initial values for the scoring algorithm
+Description:        
+        Fit an ARCH regression model to the time series y using the scoring algorithm in Engle’s original ARCH paper.
+        The model is
+        y(t) = b(1) * x(t,1) + … + b(k) * x(t,k) + e(t),
+        h(t) = a(1) + a(2) * e(t-1)^2 + … + a(p+1) * e(t-p)^2
+        in which e(t) is N(0, h(t)), given a time-series vector y up to time t-1 and a matrix of (ordinary) regressors x upto t. The order of the regression of the residual variance is specified by p.
+        If invoked as arch_fit (y, k, p) with a positive integer k, fit an ARCH(k, p) process, i.e., do the above with the t-th row of x given by
+        [1, y(t-1), …, y(t-k)]
+        Optionally, one can specify the number of iterations iter, the updating factor gamma, and initial values a0 and b0 for the scoring algorithm.
+*/
+function [a, b] = arch_fit (y, x, p, iter, gamma, a0, b0)
+    nargin = argn(2)
+    if (nargin < 3 || nargin == 6)
+      error("invalid inputs");
+    end
+    if (~ (isvector (y)))
+      error ("arch_fit: Y must be a vector");
+    end
+    T = max(size(y));
+    y = matrix (y, T, 1);
+    [rx, cx] = size (x);
+    if ((rx == 1) && (cx == 1))
+      x = autoreg_matrix (y, x);
+    elseif (~ (rx == T))
+      error ("arch_fit: either rows (X) == length (Y), or X is a scalar");
+    end
+    [T, k] = size (x);
+    if (nargin == 7)
+      a = a0;
+      b = b0;
+      e = y - x * b;
+    else
+      [b, v_b, e] = ols (y, x);
+      zer = zeros(1,p);
+      a = [v_b zer]';
+      if (nargin < 5)
+        gamma = 0.1;
+        if (nargin < 4)
+          iter = 50;
+        end
+      end
+    end
+    esq = e.^2;
+    Z = autoreg_matrix (esq, p);
+    for i = 1 : iter
+      h   = Z * a;
+      tmp = esq ./ h.^2 - 1 ./ h;
+      s   = 1 ./ h(1:T-p);
+      for j = 1 : p
+        s = s - a(j+1) * tmp(j+1:T-p+j);
+      end
+      r = 1 ./ h(1:T-p);
+      for j = 1:p
+        r = r + 2 * h(j+1:T-p+j).^2 .* esq(1:T-p);
+      end
+      r = sqrt (r);
+      X_tilde = x(1:T-p, :) .* (r * ones (1,k));
+      e_tilde = e(1:T-p) .*s ./ r;
+      delta_b = inv (X_tilde' * X_tilde) * X_tilde' * e_tilde;
+      b  = b + gamma * delta_b;  
+      e   = y - x * b;
+      esq = e .^ 2;
+      if isempty(esq) then
+          esq = zeros(size(y))
+      end
+      Z   = autoreg_matrix (esq, p);
+      h   = Z * a;
+      f   = esq ./ h - ones (T,1);
+      Z_tilde = Z ./ (h * ones (1, p+1));
+      delta_a = inv (Z_tilde' * Z_tilde) * Z_tilde' * f;
+      a = a + gamma * delta_a;
+    end
+  endfunction