close Warning: Can't synchronize with repository "(default)" (/var/svn/mms does not appear to be a Subversion repository.). Look in the Trac log for more information.

Ticket #904: _general.tol

File _general.tol, 7.3 KB (added by lmperez, 14 years ago)
Line 
1//////////////////////////////////////////////////////////////////////////////
2// FILE    : _general.tol
3// AUTHOR : Luis Miguel Pérez Liberal (lmperez@bayesinf.com)
4// PURPOSE : Standard functions for diagnosis procedures
5//////////////////////////////////////////////////////////////////////////////
6
7
8//////////////////////////////////////////////////////////////////////////////
9Set DatesIndComplet(Serie ind)
10//////////////////////////////////////////////////////////////////////////////
11{
12  Set Select(Dates(Dating(ind),First(ind),Last(ind)),Real (Date dte){
13    Real If(BinEQ(SerDat(ind,dte),?),False,True)
14  })
15};
16//////////////////////////////////////////////////////////////////////////////
17PutDescription("Dada una serie, devuelve el conjunto de días donde está
18definida por su fechado, aunque éstos sean cero",
19DatesIndComplet);
20//////////////////////////////////////////////////////////////////////////////
21
22
23//////////////////////////////////////////////////////////////////////////////
24Polyn GetRandomPolyn(Polyn p)
25//////////////////////////////////////////////////////////////////////////////
26{
27  Polyn SetSum(EvalSet(Monomes(p),Polyn (Polyn m)
28  {
29    Polyn If(EQ(Coef(m,0),1), Polyn 1, Polyn Rand(0.05,0.95)*B^Degree(m))
30  }))
31};
32//////////////////////////////////////////////////////////////////////////////
33PutDescription("Generates a random polyn with values between 0.05 and 0.95",
34GetRandomPolyn);
35//////////////////////////////////////////////////////////////////////////////
36
37
38//////////////////////////////////////////////////////////////////////////////
39Real CalcArimaParms(Text arimaLabel)
40//////////////////////////////////////////////////////////////////////////////
41{
42  Set sArima  = GetArimaFromLabel(arimaLabel);
43  Polyn p.AR  = SetProd(Traspose(sArima)[2]);
44  Polyn p.MA  = SetProd(Traspose(sArima)[3]);
45
46  Real SetSum(EvalSet(Monomes(p.AR),Real (Polyn m){If(EQ(Coef(m,0),1),0,1)}))+
47  Real SetSum(EvalSet(Monomes(p.MA),Real (Polyn m){If(EQ(Coef(m,0),1),0,1)}))
48};
49//////////////////////////////////////////////////////////////////////////////
50PutDescription("Calcula el numero de parametros arima de una arima Label",
51CalcArimaParms);
52//////////////////////////////////////////////////////////////////////////////
53
54
55//////////////////////////////////////////////////////////////////////////////
56Serie GenSerie.Normal(Anything x, Text arimaLabel)
57//////////////////////////////////////////////////////////////////////////////
58{
59  Set sArima  = GetArimaFromLabel(arimaLabel);
60  Polyn p.AR  = SetProd(Traspose(sArima)[2]);
61  Polyn p.MA  = SetProd(Traspose(sArima)[3]);
62  Polyn p.DIF = SetProd(Traspose(sArima)[4]);
63  Polyn ar    = GetRandomPolyn(p.AR);
64  Polyn ma    = GetRandomPolyn(p.MA);
65  Ratio psi   = ma/(p.DIF*ar);
66
67  Serie at = If(Grammar(x)=="Serie", x,{
68    Serie SubSer(Gaussian(0,Max(0.5,Gaussian(1,2))),Today-x/2,Today+x/2)
69  });
70
71  Serie Exp(Log(DifEq(psi,at,100)))
72};
73//////////////////////////////////////////////////////////////////////////////
74PutDescription("Generates an aleatory series without filter.
75Receives a residuals serie or a Real with longitude of days.
76Receives also an arima label to generate a time structure.
77
78Examples:
79Serie GenSerie.Normal(GenResiduals.Normal(500),\"P1DIF1AR0MA1\");
80Serie GenSerie.Normal(500,\"P1DIF1AR0MA1\");
81",
82GenSerie.Normal);
83//////////////////////////////////////////////////////////////////////////////
84
85
86//////////////////////////////////////////////////////////////////////////////
87Serie GenSerie.Heteroscedastic(Anything x, Text arimaLabel)
88//////////////////////////////////////////////////////////////////////////////
89{
90  Set sArima  = GetArimaFromLabel(arimaLabel);
91  Polyn p.AR  = SetProd(Traspose(sArima)[2]);
92  Polyn p.MA  = SetProd(Traspose(sArima)[3]);
93  Polyn p.DIF = SetProd(Traspose(sArima)[4]);
94  Polyn ar    = GetRandomPolyn(p.AR);
95  Polyn ma    = GetRandomPolyn(p.MA);
96  Ratio psi   = ma/(p.DIF*ar);
97
98  Serie at = If(Grammar(x)=="Serie", x,{
99    Real  day  = Round(Rand(x/2*0.1,x/2*0.9));
100    Serie at1  = SubSer(Gaussian(0,Max(0.5,Gaussian(1,2))),Today-x/2,Today-day);
101    Serie at2  = SubSer(Gaussian(0,Max(0.5,Gaussian(1,2))),Today-day, Today+x/2);
102    Serie SigK = SetSumC(For(1,7,Serie (Real j){
103      Serie CalInd(WD(j),Dating(at1))*Max(0.5,Gaussian(1,2))
104    }));
105    Serie SetSumC([[ at1,at2 ]])*SigK
106  });
107
108  Serie Exp(Log(DifEq(psi,at,100)))
109};
110//////////////////////////////////////////////////////////////////////////////
111PutDescription("Generates an aleatory series with heteroscedasticity and with
112'long' days. Receives an arima label to generate a time structure.
113
114Examples:
115Serie GenSerie.Heteroscedastic(GenResiduals.Heteroscedastic(500),\"P1DIF1AR0MA1\");
116Serie GenSerie.Normal(GenResiduals.Heteroscedastic(500),\"P1DIF1AR0MA1\");
117Serie GenSerie.Heteroscedastic(500,\"P1DIF1AR0MA1\");",
118GenSerie.Heteroscedastic);
119//////////////////////////////////////////////////////////////////////////////
120
121
122//////////////////////////////////////////////////////////////////////////////
123Serie GenResiduals.Normal(Real long)
124//////////////////////////////////////////////////////////////////////////////
125{
126  Serie SubSer(Gaussian(0,Max(0.5,Gaussian(1,2))),Today-long/2,Today+long/2)
127};
128//////////////////////////////////////////////////////////////////////////////
129PutDescription("Generates an aleatory residuals serie with 'long' days",
130GenResiduals.Normal);
131//////////////////////////////////////////////////////////////////////////////
132
133
134//////////////////////////////////////////////////////////////////////////////
135Serie GenResiduals.Heteroscedastic(Real long)
136//////////////////////////////////////////////////////////////////////////////
137{
138  Real  day = Round(Rand(long/2*0.2,long/2*0.8));
139  Serie at1 = SubSer(Gaussian(0,Max(0.5,Gaussian(1,2))),Today-long/2,Today-day);
140  Serie at2 = SubSer(Gaussian(0,Max(0.5,Gaussian(1,2))),Today-day, Today+long/2);
141  Serie SigK = SetSumC(For(1,7,Serie (Real j){
142    Serie CalInd(WD(j),Dating(at1))*Max(0.5,Gaussian(1,2))
143  }));
144
145  Serie SetSumC([[ at1,at2 ]])*SigK
146};
147//////////////////////////////////////////////////////////////////////////////
148PutDescription("Generates an aleatory residuals serie with heteroscedasticity
149and with 'long' days",
150GenResiduals.Heteroscedastic);
151//////////////////////////////////////////////////////////////////////////////
152
153
154//////////////////////////////////////////////////////////////////////////////
155Real BinEQ(Real x, Real y)
156//////////////////////////////////////////////////////////////////////////////
157{
158  Real Case( And(     IsUnknown(x),      IsUnknown(y)  ),  Real True,
159             And(     IsUnknown(x),  Not(IsUnknown(y)) ),  Real False,
160             And( Not(IsUnknown(x)),     IsUnknown(y)  ),  Real False,
161             And( Not(IsUnknown(x)), Not(IsUnknown(y)) ),  Real EQ(x,y))
162};
163//////////////////////////////////////////////////////////////////////////////
164PutDescription("Binary operator like EQ or ==, between real numbers,
165but considering a ommit value like number (null), not like an irresolution",
166BinEQ);
167//////////////////////////////////////////////////////////////////////////////
168
169
170//////////////////////////////////////////////////////////////////////////////
171// END
172//////////////////////////////////////////////////////////////////////////////