Ticket #1608: ejemplo1.tol

//////////////////////////////////////////////////////////////////////////////
// File: ejemplo1.tol
// Purpuse: ejemplo inventado de un modelo mixto
//////////////////////////////////////////////////////////////////////////////

//////////////////////////////////////////////////////////////////////////////
#Require MMS;

// Real PutRandomSeed(190463);
// Real PutRandomSeed(514809833);
Real PutRandomSeed(1528068215); 


// Definimos las series de las observaciones y las series que entraran
// tanto como inputs aditivos como multiplicativos.
  
Date begin = y2000m01;
Date end   = y2012m08;


Polyn ma = 1-0.3*B^12;
Polyn ar = 1-0.6*B;
Polyn dif = (1-B^12);
Real initLength = 100*Max(Degree(ar)+ Degree(dif), Degree(ma));
Date begin_ = Succ(begin, Monthly, -initLength);  
Serie R0 = SubSer(Gaussian(0, 0.2, Monthly), begin_, end); 
Serie N0 = DifEq(ma/(dif*ar), R0);
Serie R = SubSer(R0, begin, end);
Serie N = SubSer(N0, begin, end)+4;
  
// X1 ~ N(10, 2)
Real beta1 = 2;
Serie X1 = SubSer(Gaussian(3, 1, Monthly), begin, end);

// X2 ~ N(5, 4)
Real beta2 = 0.2;
Serie X2_ = SubSer(Gaussian(5, 0.4, Monthly), begin, end);
Serie X2 = DifEq((1)/(1-0.3*B), X2_, N0);

// X3 ~ Festivos
TimeSet CtFes.es = 
  M(1)*D(1) + // Año nuevo
  M(1)*D(6) + // Epifania (no siempre es FasNac)
  Succ(Easter, -2, C) + // Viernes Santo
  M(5)*D(1)   + // Dia del trabajador
  M(8)*D(15)  + // Virgen de Agosto
  M(10)*D(12) + // Dia del Pilar
  M(11)*D(1)  + // Todos los Santos
  M(12)*D(6)  + // Constitucion
  M(12)*D(8)  + // Inmaculada
  M(12)*D(25) +  // Navidad
  WD(7) + WD(6)
;

Real beta3 = 0.20;
Serie X3 = SubSer(CalInd(CtFes.es, Monthly), begin, end);
  
Serie additive = beta1*X1;
Serie multiplicative = N + beta2*X2 + beta3*X3; 
Serie observation = Exp(multiplicative)+additive;

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



//////////////////////////////////////////////////////////////////////////////
// DATASET: Conjunto de variables

MMS::@DataSet DS = MMS::Container::ReplaceDataSet([[
  Text _.name = "Mixto_DataSet"
]]);

// ReplaceVariable

Anything DS::ReplaceVariable([[
  Text _.name = "Observation";
  Text _.description = "Variable de las observaciones";
  Text _.expression = "Serie observation";
  Set _.tags = SetOfText ("Output")
]]);

Anything DS::ReplaceVariable([[
  Text _.name = "Input1";
  Text _.description = "Variable Input1";
  Text _.expression = "Serie X1";
  Set _.tags = SetOfText ("Input", "additive") 
]]);

Anything DS::ReplaceVariable([[
  Text _.name = "Input2";
  Text _.description = "Variable Input2";
  Text _.expression = "Serie X2";
  Set _.tags = SetOfText ("Input", "multiplicative")
]]);

Anything DS::ReplaceVariable([[
  Text _.name = "Input3";
  Text _.description = "Variable Input3";
  Text _.expression = "Serie X3";
  Set _.tags = SetOfText ("Input", "multiplicative")
]]);
//////////////////////////////////////////////////////////////////////////////



//////////////////////////////////////////////////////////////////////////////
// MODEL: Modelo y Submodelo: TerminosExplicativos+ Noise

// CreateModel
MMS::@Model model = MMS::Container::ReplaceModel([[
  Text _.name = "Mixto_Model";
  Text _.description = "Modelo mixto: aditivo y multiplicativo";
  Set _.dataSets = [[ "Mixto_DataSet" ]]
]]);

// CreateSubmodel
MMS::@Submodel submodel = model::CreateSubmodel([[
  Text _.name = "Mixto_Submodel";
  Text _.descripton = "Submodelo del ejemplo de modelo mixto";
  NameBlock _.output = [[
    Text _.name = "Output";
    Text _.variable = "Observation";
    Text _.transformationLabel = "BoxCox_0_0"
  ]];
  NameBlock _.noise = [[
    Text _.type = "ARIMA"; 
    Text _.arimaLabel = "P1_12DIF0_1AR1_0MA0_12";
    Real _.sigma = 1
  ]]
]]);

// Para crear el input aditivo, simplemente tenemos que poner el 
// atributo de la clase término explicativo:
//
// Real _.isAdditive = 1 (True)
//
// Por defecto, este atributo está con valor 0 (False)
// Observacion
// el polinomio de la función de transferencia no podemos introducir un 0
// porque no lo toma como polinomio constante 0, lo toma como que no 
// existe polinomio.

// CreateExpTerm_TransferFunction
Set EvalSet(Select(
  model::GetDataSet(?)::GetVariables(?), Real(NameBlock var){
    "Input" <: var::GetTags(?)
  }), 
  Anything (NameBlock varInput){
    submodel::CreateExpTerm_TransferFunction([[
      Text _.name  = varInput::GetName(?);
      Real _.isAdditive = If("additive" <: varInput::GetTags(?), 1, 0);
      NameBlock _.input = [[
        Text _.name = varInput::GetName(?);
        Text _.variable = varInput::GetName(?)
      ]];
      Polyn _.transferFunction = 0.1
    ]])
});

// Restricción del término aditivo
/*
Real model::GetParameter("Mixto_Submodel__Input1__Linear.0")
  ::SetPrior([[
    Real _.mean = 2;
    Real _.sigma = 0.1
]]);
*/
//////////////////////////////////////////////////////////////////////////////



//////////////////////////////////////////////////////////////////////////////
// ESTIMATION: Estimación MLE
MMS::@Estimation estimationMLE= MMS::Container::ReplaceEstimation([[
  Text _.name = "Mixto_Estimate";
  MMS::@Model _.model = MMS::Container::GetModel("Mixto_Model");
  MMS::@SettingsMultiMLE _.settings = [[
    Real _.showTraces = False  // Para que no muestre la traza
  ]]
]]);

// Ejecución de la estimación:
Real estimationMLE::Execute(?);


// ESTIMATION: Estimación BSR
MMS::@Estimation estimationBSR = MMS::Container::ReplaceEstimation([[
  Text _.name = "Mixto_EstimateBSR";
  MMS::@Model _.model = MMS::Container::GetModel("Mixto_Model");
  MMS::@SettingsBSR _.settings = [[
    Real mcmc.sampleLength = 500; // tamaño de la muestra
    Real do.report = True 
  ]]
]]);

// Ejecución de la estimación:
Real estimationBSR::Execute(?);

// ESTIMATION: Estimación BSR
MMS::@Estimation estimationBSR2 = MMS::Container::ReplaceEstimation([[
  Text _.name = "Mixto_EstimateBSR2";
  MMS::@Model _.model = MMS::Container::GetModel("Mixto_Model");
  MMS::@SettingsBSR _.settings = [[
    Real _.additiveApproximation = 0;
    Real mcmc.sampleLength = 500; // tamaño de la muestra
    Real do.report = True 
  ]]
]]);

// Ejecución de la estimación:
Real estimationBSR2::Execute(?);
//////////////////////////////////////////////////////////////////////////////