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 #167: ej_cov_3.tol

File ej_cov_3.tol, 4.2 KB (added by Pedro Gea, 13 years ago)
Line 
1#Require MMS;
2Real PutRandomSeed(2143);
3
4Date begin0 = y1989m12;
5Date begin = y2000;
6Date end = y2012m12;
7Real phi = 0.5;
8Serie residuals = SubSer(Gaussian(0, 0.5, Monthly), begin0, end);
9Serie noise = SubSer(DifEq(1/(1-phi*B), residuals, 0), begin, end);
10Serie input1_ = SubSer(Rand(-1, 1, Monthly), begin, end);
11Serie input2_ = SubSer(Floor(Rand(0, 2, Monthly)), begin, end);
12Real beta1 = 0.4;
13Real beta2 = -0.2;
14Serie output_ = noise + beta1*input1_ + beta2*input2_;
15
16MMS::@DataSet dataSet = MMS::Container::ReplaceDataSet([[
17  Text _.name = "DataSet"
18]]);
19Anything dataSet::CreateVariable_Data(Serie {
20  Serie output = Copy(output_)
21});
22Anything dataSet::CreateVariable_Data(Serie {
23  Serie input1 = Copy(input1_)
24});
25Anything dataSet::CreateVariable_Data(Serie {
26  Serie input2 = Copy(input2_)
27});
28
29//////////////////////////////////////////////////////////////////////////////
30
31MMS::@Model model = MMS::Container::ReplaceModel([[
32  Text _.name = "Model1";
33  MMS::@DataSet _.dataSet = dataSet
34]]);
35
36MMS::@Submodel submodel = model::CreateSubmodel([[
37  Text _.name = "Submodel";
38  NameBlock _.output = [[
39    Text _.name = "output";
40    Text _.variable = "output"
41  ]];
42  NameBlock _.noise =  [[
43    Text _.type = "Normal";
44    Real _.sigma = 0.5;
45    Real _.sigmaFixed = 1
46  ]]
47]]);
48
49Anything submodel::CreateExpTerm_Coefficient([[
50  Text _.name = "ExpTerm1";
51  NameBlock _.input = [[
52    Text _.name = "input1";
53    Text _.variable = "input1"
54  ]];
55  Real _.coefficient = 0.1
56]]);
57
58Anything submodel::CreateExpTerm_Coefficient([[
59  Text _.name = "ExpTerm2";
60  NameBlock _.input = [[
61    Text _.name = "input2";
62    Text _.variable = "input2"
63  ]];
64  Real _.coefficient = -0.1
65]]);
66
67MMS::@Estimation estimation = MMS::Container::ReplaceEstimation([[
68  Text _.name = "Estimation1a";
69  MMS::@Model _.model = model;
70  MMS::@SettingsBSR _.settings = [[
71    Real _.showTraces = False;
72    Real mcmc.burnin = 100;
73    Real mcmc.sampleLength = 500
74  ]]
75]]);
76
77//////////////////////////////////////////////////////////////////////////////
78
79MMS::@Model model2 = model::Duplicate(?);
80
81Real model2::GetSubmodel(1)::SetNoise([[
82  Text _.type = "ARIMA";
83  Text _.arimaLabel = "P1DIF0AR1MA0";
84  Real _.sigma = 0.5;
85  Real _.sigmaFixed = 1
86]]);
87
88MMS::@Noise noise2 = model2::GetSubmodel(1)::GetNoise(?);
89
90MMS::@Parameter ar = noise2::GetARIMABlock(1)::GetParameterARIMA(1);
91Real ar::SetInitialValue(0.5);
92Real ar::SetIsFixed(1);
93
94MMS::@Estimation estimation2 = MMS::Container::ReplaceEstimation([[
95  Text _.name = "Estimation1b";
96  MMS::@Model _.model = model2;
97  MMS::@SettingsBSR _.settings = [[
98    Real _.showTraces = False;
99    Real mcmc.burnin = 100;
100    Real mcmc.sampleLength = 500
101  ]]
102]]);
103
104//////////////////////////////////////////////////////////////////////////////
105
106MMS::@Model model3 =  model::Duplicate(?);
107
108Real size = CountS(output_);
109Polyn pAr = Expand(1/(1-phi*B), size+100);
110Matrix cov1 = Tra(PolMat(pAr,size,size)) * PolMat(pAr,size,size);
111Matrix cov2 = PolMat(pAr,size,size) * Tra(PolMat(pAr,size,size));
112Matrix cov3 = Tra(PolMat(pAr,size+100,size)) * PolMat(pAr,size+100,size);
113
114Real model3::GetSubmodel(1)::SetNoise([[
115  Text _.type = "Normal";
116  Matrix _.relativeCovariance = cov3;
117  Real _.sigma = 0.5;
118  Real _.sigmaFixed = 1
119]]);
120
121MMS::@Estimation estimation3 = MMS::Container::ReplaceEstimation([[
122  Text _.name = "Estimation1c";
123  MMS::@Model _.model = model3;
124  MMS::@SettingsBSR _.settings = [[
125    Real _.showTraces = False;
126    Real mcmc.burnin = 100;
127    Real mcmc.sampleLength = 500
128  ]]
129]]);
130
131//////////////////////////////////////////////////////////////////////////////
132
133Real ShowResults(MMS::@Estimation est) {
134  Set EvalSet(est::GetParameters(?), Real (Anything r) {
135    If(Grammar(r)=="Real",
136      WriteLn(Name(r)<<" "<<r),
137      WriteLn(Name(r)<<" "<<r::GetMean(?)) );
138  1});
1391};
140
141Real PutRandomSeed(Time);
142
143Real PutRandomSeed(2143);
144Real estimation::Execute(?);
145Real PutRandomSeed(2143);
146Real estimation2::Execute(?);
147Real PutRandomSeed(2143);
148Real estimation3::Execute(?);
149
150Real ShowResults(estimation);
151Real ShowResults(estimation2);
152Real ShowResults(estimation3);
153
154