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
|
//Return a simulation of the ARMA model.
//Calling Sequence
//arma_rnd (a, b, v, t, n)
//arma_rnd (a, b, v, t)
//Parameters
//a: vector
//b: vector
//v: Variance
//t: Length of output vector
//n: Number of dummy x(i) used for initialization
//Description
//This is an Octave function.
//The ARMA model is defined by
//
//x(n) = a(1) * x(n-1) + … + a(k) * x(n-k)
// + e(n) + b(1) * e(n-1) + … + b(l) * e(n-l)
//in which k is the length of vector a, l is the length of vector b and e is Gaussian white noise with variance v. The function returns a vector of length t.
//
//The optional parameter n gives the number of dummy x(i) used for initialization, i.e., a sequence of length t+n is generated and x(n+1:t+n) is returned. If n is omitted, n = 100 is used.
//Examples
//a = [1 2 3 4 5];
//b = [7 8 9 10 11];
//v = 10;
//t = 5;
//n = 100;
//arma_rnd (a, b, v, t, n)
//Output :
// ans =
//
// 61400.907
// 158177.11
// 407440.29
// 1049604.
// 2703841.3
//function res = arma_rnd (a, b, v, t, n)
//funcprot(0);
//lhs = argn(1)
//rhs = argn(2)
//if (rhs < 5 | rhs > 6)
//error("Wrong number of input arguments.")
//end
//
//select(rhs)
//
// case 5 then
// res = callOctave("arma_rnd",a, b, v, t)
//
// case 6 then
// res = callOctave("arma_rnd",a, b, v, t, n)
//
// end
//endfunction
function x = arma_rnd (a, b, v, t, n)
funcprot();
[nargout,nargin] = argn() ;
if (nargin == 4)
n = 100;
elseif (nargin == 5)
if (~ isscalar (n))
error ("arma_rnd: N must be a scalar");
end
else
error("arma_rnd: invalid input");
end
if ((min (size (a)) > 1) | (min (size (b)) > 1))
error ("arma_rnd: A and B must not be matrices");
end
if (~ isscalar (t))
error ("arma_rnd: T must be a scalar");
end
ar = length (a);
br = length (b);
a = matrix (a, ar, 1);
b = matrix (b, br, 1);
// Apply our notational convention.
a = [1; -a];
b = [1; b];
n = min (n, ar + br);
e = sqrt (v) * rand(t + n, 1);
x = filter (b, a, e);
x = x(n + 1 : t + n);
endfunction
|
