Ticket #167: ej_cov.tol

#Require MMS;
Real PutRandomSeed(2143);

Date begin = y2000;
Date end = y2012m12;
Serie noise = SubSer(Gaussian(0, 0.5, Monthly), begin, end);
Serie input1_ = SubSer(Rand(-1, 1, Monthly), begin, end);
Serie input2_ = SubSer(Floor(Rand(0, 2, Monthly)), begin, end);
Real beta1 = 0.4;
Real beta2 = -0.2;
Serie output_ = noise + beta1*input1_ + beta2*input2_;

Matrix m = ((1,1), (2, -1));
Matrix mn = Col(0, 1);
Matrix mmI =  GaussInverse(Tra(m)*m);
Real ms = 0.1;
Matrix mn2 = mmI * Tra(m) * mn;

MMS::@DataSet dataSet = MMS::Container::ReplaceDataSet([[
  Text _.name = "DataSet"
]]);
Anything dataSet::CreateVariable_Data(Serie {
  Serie output = Copy(output_)
});
Anything dataSet::CreateVariable_Data(Serie {
  Serie input1 = Copy(input1_)
});
Anything dataSet::CreateVariable_Data(Serie {
  Serie input2 = Copy(input2_)
});

//////////////////////////////////////////////////////////////////////////////

MMS::@Model model = MMS::Container::ReplaceModel([[
  Text _.name = "Model1";
  MMS::@DataSet _.dataSet = dataSet
]]);

MMS::@Submodel submodel = model::CreateSubmodel([[
  Text _.name = "Submodel";
  NameBlock _.output = [[
    Text _.name = "output";
    Text _.variable = "output"
  ]];
  NameBlock _.noise =  [[
    Text _.type = "Normal"
  ]]
]]);
   
Anything submodel::CreateExpTerm_Coefficient([[
  Text _.name = "ExpTerm1";
  NameBlock _.input = [[
    Text _.name = "input1";
    Text _.variable = "input1"
  ]];
  Real _.coefficient = 0.1
]]);

Anything submodel::CreateExpTerm_Coefficient([[
  Text _.name = "ExpTerm2";
  NameBlock _.input = [[
    Text _.name = "input2";
    Text _.variable = "input2"
  ]];
  Real _.coefficient = -0.1
]]);

MMS::@Estimation estimation = MMS::Container::ReplaceEstimation([[
  Text _.name = "Estimation1a";
  MMS::@Model _.model = model;
  MMS::@SettingsBSR _.settings = [[
    Real _.showTraces = False;
    Real mcmc.burnin = 100;
    Real mcmc.sampleLength = 500
  ]]
]]);

//////////////////////////////////////////////////////////////////////////////

MMS::@Model model2 =  model::Duplicate(?);

Anything model2::CreateMCombination([[
  Text _.name = "Suma";
  Set _.parameters = [[
    "Submodel__ExpTerm1__Linear.0";
    "Submodel__ExpTerm2__Linear.0"
  ]];
  Set _.coefficients = [[MatDat(m,1,1), MatDat(m,1,2)]];
  NameBlock _.prior = [[
    Real _.mean = MatDat(mn,1,1);
    Real _.sigma = ms
  ]]
]]);
Anything model2::CreateMCombination([[
  Text _.name = "Diferencia";
  Set _.parameters = [[
    "Submodel__ExpTerm1__Linear.0";
    "Submodel__ExpTerm2__Linear.0"
  ]];
  Set _.coefficients = [[MatDat(m,2,1), MatDat(m,2,2)]];
  NameBlock _.prior = [[
    Real _.mean = MatDat(mn,2,1);
    Real _.sigma = ms
  ]]
]]);

MMS::@Estimation estimation2 = MMS::Container::ReplaceEstimation([[
  Text _.name = "Estimation1b";
  MMS::@Model _.model = model2;
  MMS::@SettingsBSR _.settings = [[
    Real _.showTraces = False;
    Real mcmc.burnin = 100;
    Real mcmc.sampleLength = 500
  ]]
]]);

//////////////////////////////////////////////////////////////////////////////

MMS::@Model model3 =  model::Duplicate(?);

MMS::@Hierarchy hierarchy = model3::CreateHierarchy([[
  Text _.name = "Hierarchy";
  Set _.mElements = [[
    "Submodel__ExpTerm1__Linear.0";
    "Submodel__ExpTerm2__Linear.0"
  ]];
  NameBlock _.noise = [[
    Text _.type = "Normal";
    Real _.sigma = ms;
    Real _.sigmaFixed = 1;
    Matrix _.relativeCovariance = mmI
  ]]
]]);
Anything hierarchy::CreateHierarchyTerm([[
  Text _.name = "Term";
  Real _.initialParameter = 1;
  Real _.isFixed = 1;
  Set _.coefficients = [[ MatDat(mn2,1,1), MatDat(mn2,2,1) ]]  
]]);

MMS::@Estimation estimation3 = MMS::Container::ReplaceEstimation([[
  Text _.name = "Estimation1c";
  MMS::@Model _.model = model3;
  MMS::@SettingsBSR _.settings = [[
    Real _.showTraces = False;
    Real mcmc.burnin = 100;
    Real mcmc.sampleLength = 500
  ]]
]]);

//////////////////////////////////////////////////////////////////////////////

MMS::@Model model4 =  model::Duplicate(?);

Anything model4::CreateMCombination([[
  Text _.name = "Suma";
  Set _.parameters = [[
    "Submodel__ExpTerm1__Linear.0";
    "Submodel__ExpTerm2__Linear.0"
  ]];
  Set _.coefficients = [[MatDat(m,1,1), MatDat(m,1,2)]]
]]);
Anything model4::CreateMCombination([[
  Text _.name = "Diferencia";
  Set _.parameters = [[
    "Submodel__ExpTerm1__Linear.0";
    "Submodel__ExpTerm2__Linear.0"
  ]];
  Set _.coefficients = [[MatDat(m,2,1), MatDat(m,2,2)]]
]]);

MMS::@Hierarchy hierarchy2 = model4::CreateHierarchy([[
  Text _.name = "Hierarchy";
  Set _.mElements = [[
    "Suma";
    "Diferencia"
  ]];
  NameBlock _.noise = [[
    Text _.type = "Normal";
    Real _.sigma = ms;
    Real _.sigmaFixed = 1
  ]]
]]);
Anything hierarchy2::CreateHierarchyTerm([[
  Text _.name = "Term";
  Real _.initialParameter = 1;
  Real _.isFixed = 1;
  Set _.coefficients = [[ MatDat(mn,1,1), MatDat(mn,2,1) ]]
]]);

MMS::@Estimation estimation4 = MMS::Container::ReplaceEstimation([[
  Text _.name = "Estimation1d";
  MMS::@Model _.model = model4;
  MMS::@SettingsBSR _.settings = [[
    Real _.showTraces = False;
    Real mcmc.burnin = 100;
    Real mcmc.sampleLength = 500
  ]]
]]);

//////////////////////////////////////////////////////////////////////////////

Real ShowResults(MMS::@Estimation est) {
  Set EvalSet(est::GetParameters(?), Real (Anything r) {
    If(Grammar(r)=="Real",
      WriteLn(Name(r)<<" "<<r),
      WriteLn(Name(r)<<" "<<r::GetMean(?)) );
  1});
1};


Real PutRandomSeed(Time);

Real PutRandomSeed(2143);
Real estimation::Execute(?);
Real PutRandomSeed(2143);
Real estimation2::Execute(?);
Real PutRandomSeed(2143);
Real estimation3::Execute(?);
Real PutRandomSeed(2143);
Real estimation4::Execute(?);

Real ShowResults(estimation);
Real ShowResults(estimation2);
Real ShowResults(estimation3);
Real ShowResults(estimation4);