It is some time ago, that I was fighting with that model and how
to estimate the 3 parameters (until recognizing that one may want
restrictions for them).
Essentially one uses the statistics of the underlying data to get
a reasonable starting guess. 
For an estimation then one can use the Optimization package.
The gradients involved are best coded as a floating point library,
which here is done through a DLL (it should work on all Windows OS),
code in C is included for that.
Meanwhile all can be done using 'Compiler:-compile' or the Watcom
compiler delivered with Maple.
I left that older sheet as it is - in Maple 9 (may be one should
brush it up for concurrent Maple versions, they are grumbling a bit
about the code) and hope it is of interest despite of that.

Maple sheet: www.mapleprimes.com/files/102_Garch.mws
DLL + code:  www.mapleprimes.com/files/102_Garch.dll.zip
as a pdf:    www.mapleprimes.com/files/102_Garch.mws.pdf

  Int(phi(u-a*I)*exp(k*u*I),u = -infinity .. infinity);
  Change(%, u-a*I = s, s); subs(s=u,%): expand(%): combine(%,exp);

                   infinity
                  /
                 |
   ...

  restart; interface(version);

    Classic Worksheet Interface, Maple 12.00, Windows, Apr 10 2008 Build ID 347164

  Int(1/(u^2+1/4)*exp(I*u*k),u=-infinity...infinity); value(%);

                                  0

but

  Int(1/(u^2+1/4)*exp(I*u*k),u=-infinity...infinity...

This behaviour of 'assume' puzzles me ... 
and it needed some time to locate it
 
  restart; interface(version);

  Classic Worksheet Interface, Maple 12.00, Windows, Apr 10 2008 Build ID 347164

  assume(a::real); 
    #a:='a'; # activating removes the problem
  f0:=x+a; 
    #a:='a'; # that would not help - since is in use?
  assume(0<a);
  f:=collect(f0,x);

     ...

Here is some standard alternative to Newton's method (and
thus may safe some homework ... so what). 

It will find a root of f (I think it must be continuous C^1)
in the interval ax ... bx, if it has different signs at the
boundaries.

The code is more or less translated from netlib C library
(or similar).

Usage:

  Digits:=16;

  f:= x -> exp(x)-Pi;
  zBrent( f, 0.0, 2.0);
  'f(%)': '%'= evalf(%);