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_2.tol

File ej_cov_2.tol, 4.0 KB (added by Pedro Gea, 13 years ago)
Line 
1#Require MMS;
2Real PutRandomSeed(2143);
3
4Date begin0 = y1999m12;
5Date begin = y2000;
6Date end = y2012m12;
7Real phi = 0.5;
8Serie residuals = SubSer(Gaussian(0, 0.5, Monthly), begin0, end);
9Serie noise = (1-phi*B):residuals;
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  ]]
45]]);
46
47Anything submodel::CreateExpTerm_Coefficient([[
48  Text _.name = "ExpTerm1";
49  NameBlock _.input = [[
50    Text _.name = "input1";
51    Text _.variable = "input1"
52  ]];
53  Real _.coefficient = 0.1
54]]);
55
56Anything submodel::CreateExpTerm_Coefficient([[
57  Text _.name = "ExpTerm2";
58  NameBlock _.input = [[
59    Text _.name = "input2";
60    Text _.variable = "input2"
61  ]];
62  Real _.coefficient = -0.1
63]]);
64
65MMS::@Estimation estimation = MMS::Container::ReplaceEstimation([[
66  Text _.name = "Estimation1a";
67  MMS::@Model _.model = model;
68  MMS::@SettingsBSR _.settings = [[
69    Real _.showTraces = False;
70    Real mcmc.burnin = 100;
71    Real mcmc.sampleLength = 500
72  ]]
73]]);
74
75
76//////////////////////////////////////////////////////////////////////////////
77
78MMS::@Model model2 = model::Duplicate(?);
79
80Real model2::GetSubmodel(1)::SetNoise([[
81  Text _.type = "ARIMA";
82  Text _.arimaLabel = "P1DIF0AR0MA1"
83]]);
84
85MMS::@Noise noise2 = model2::GetSubmodel(1)::GetNoise(?);
86
87MMS::@Parameter ma = noise2::GetARIMABlock(1)::GetParameterARIMA(1);
88Real ma::SetInitialValue(0.5);
89Real ma::SetIsFixed(1);
90
91MMS::@Estimation estimation2 = MMS::Container::ReplaceEstimation([[
92  Text _.name = "Estimation1b";
93  MMS::@Model _.model = model2;
94  MMS::@SettingsBSR _.settings = [[
95    Real _.showTraces = False;
96    Real mcmc.burnin = 100;
97    Real mcmc.sampleLength = 500
98  ]]
99]]);
100
101//////////////////////////////////////////////////////////////////////////////
102
103MMS::@Model model3 =  model::Duplicate(?);
104
105Real size = CountS(output_);
106Polyn pMa = 1-phi*B;
107Matrix cov1 = Tra(PolMat(pMa,size,size)) * PolMat(pMa,size,size);
108Matrix cov2 = PolMat(pMa,size,size) * Tra(PolMat(pMa,size,size));
109Matrix cov3 = Tra(PolMat(pMa,size+1,size)) * PolMat(pMa,size+1,size);
110
111Real model3::GetSubmodel(1)::SetNoise([[
112  Text _.type = "Normal";
113  Matrix _.relativeCovariance = cov3
114]]);
115
116MMS::@Estimation estimation3 = MMS::Container::ReplaceEstimation([[
117  Text _.name = "Estimation1c";
118  MMS::@Model _.model = model3;
119  MMS::@SettingsBSR _.settings = [[
120    Real _.showTraces = False;
121    Real mcmc.burnin = 100;
122    Real mcmc.sampleLength = 500
123  ]]
124]]);
125
126//////////////////////////////////////////////////////////////////////////////
127
128Real ShowResults(MMS::@Estimation est) {
129  Set EvalSet(est::GetParameters(?), Real (Anything r) {
130    If(Grammar(r)=="Real",
131      WriteLn(Name(r)<<" "<<r),
132      WriteLn(Name(r)<<" "<<r::GetMean(?)) );
133  1});
1341};
135
136
137Real PutRandomSeed(Time);
138
139Real PutRandomSeed(2143);
140Real estimation::Execute(?);
141Real PutRandomSeed(2143);
142Real estimation2::Execute(?);
143Real PutRandomSeed(2143);
144Real estimation3::Execute(?);
145
146Real ShowResults(estimation);
147Real ShowResults(estimation2);
148Real ShowResults(estimation3);
149