1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
|
//Simulate an ARCH sequence of length t with AR coefficients b and CH coefficients a.
//Calling Sequence
//arch_rnd (a, b, t)
//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.
//The result y(t) 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
//Examples
//a = [1 2 3 4 5];
//b = [7 8 9 10];
//t = 5 ;
//arch_rnd (a, b, t)
//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)
//
// 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
|
