Index: tol_tests/tol/Bugzilla/bug_001608/evalModel.tol =================================================================== --- tol_tests/tol/Bugzilla/bug_001608/evalModel.tol (revision b27f50548784bca303be5018af783b0f3f8a3f90) +++ tol_tests/tol/Bugzilla/bug_001608/evalModel.tol (revision b27f50548784bca303be5018af783b0f3f8a3f90) @@ -0,0 +1,21 @@ + +//Matrix param = SubRow(estim::report::mcmc,[[1]]); +Matrix param = estim::param.average; + +Real beta1_ = MatDat(param,1,6); +Real beta2_ = MatDat(param,1,1); +Real beta3_ = MatDat(param,1,2); +Real s2 = MatDat(param,1,3); +Real phi1 = MatDat(param,1,4); +Real theta12 = MatDat(param,1,5); +Serie N; +Serie N_ = Log(observation-beta1_*X1)-beta2_*X2-beta3_*X3; +Set arima = [[ + ARIMAStruct(1,1-phi1*B,1,1), + ARIMAStruct(12,1,1-theta12*B^12,dif)]]; +Matrix N__ = Tra(SerSetMat([[N_]])); +Set ev = ARIMALevinsonEval(arima, N__); + +Serie R; +Serie R_ = MatSerSet(Tra(ev::a),Dating(N),Succ(Last(N),Dating(N),1-Rows(ev::a)))[1]; + Index: tol_tests/tol/Bugzilla/bug_001608/generateModelData_0.tol =================================================================== --- tol_tests/tol/Bugzilla/bug_001608/generateModelData_0.tol (revision b27f50548784bca303be5018af783b0f3f8a3f90) +++ tol_tests/tol/Bugzilla/bug_001608/generateModelData_0.tol (revision b27f50548784bca303be5018af783b0f3f8a3f90) @@ -0,0 +1,75 @@ +////////////////////////////////////////////////////////////////////////////// +// 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); + +/* +Log(Y-beta1*X1)-(beta2*X2 + beta3*X3) = N +Log(Y-beta1*X1) = N + (beta2*X2 + beta3*X3) +Y = Exp(N + (beta2*X2 + beta3*X3))+beta1*X1 +*/ +Serie additive = beta1*X1; +Serie mult2 = beta2*X2; +Serie mult3 = beta3*X3; +Serie expMult2 = Exp(mult2); +Serie expMult3 = Exp(mult3); +Serie expNoise = Exp(N); +Serie expRes = Exp(R); +Serie multiplicative = N + beta2*X2 + beta3*X3; +Serie expMult = Exp(multiplicative); +Serie observation = expMult+additive; + + + Index: tol_tests/tol/Bugzilla/bug_001608/generateModelData_1.tol =================================================================== --- tol_tests/tol/Bugzilla/bug_001608/generateModelData_1.tol (revision b27f50548784bca303be5018af783b0f3f8a3f90) +++ tol_tests/tol/Bugzilla/bug_001608/generateModelData_1.tol (revision b27f50548784bca303be5018af783b0f3f8a3f90) @@ -0,0 +1,63 @@ + +// 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.9*B^12; +Polyn ar = 1-0.3*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.04, 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 = 5; +//Serie X1 = SubSer(Gaussian(10, 2, Monthly), begin, end); + Serie X1 = SubSer(CalVar(C,Monthly), begin, end); + +// X2 ~ N(5, 4) +Real beta2 = 0.2; +Serie X2_ = SubSer(Gaussian(0, 1, 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.05; +Serie X3 = SubSer(CalVar(CtFes.es, Monthly), begin, end); + +Serie additive = beta1*X1; +Serie mult2 = beta2*X2; +Serie mult3 = beta3*X3; +Serie expMult2 = Exp(mult2); +Serie expMult3 = Exp(mult3); +Serie expNoise = Exp(N); +Serie multiplicative = N + beta2*X2 + beta3*X3; +Serie expMult = Exp(multiplicative); +Serie observation = expMult+additive; +Serie zero = observation*0; + +////////////////////////////////////////////////////////////////////////////// + Index: tol_tests/tol/Bugzilla/bug_001608/runBsr.tol =================================================================== --- tol_tests/tol/Bugzilla/bug_001608/runBsr.tol (revision b27f50548784bca303be5018af783b0f3f8a3f90) +++ tol_tests/tol/Bugzilla/bug_001608/runBsr.tol (revision b27f50548784bca303be5018af783b0f3f8a3f90) @@ -0,0 +1,62 @@ +#Require MMS; + +//Set Include("generateModelData.tol"); +Text modelName = "estimationBSR2_0"; + +BysMcmc::@Config bsr.config = Include(modelName+"/master.bsr.config.oza")[1]; + +Real bsr.config::mcmc.burnin := 0; +Real bsr.config::mcmc.sampleLength := 500; + +Set bsr.parsed = BSR.Parse(modelName+"/master.bsr","master"); +Real ((((bsr.parsed::SubModule)[1])->LinearBlock)[1])->InitialValue := beta2; +Real ((((bsr.parsed::SubModule)[1])->LinearBlock)[2])->InitialValue := beta3; +Polyn((((((bsr.parsed::SubModule)[1])->NoiseDistrib)[1])->Arima)[1])->AR:=ar; +Polyn((((((bsr.parsed::SubModule)[1])->NoiseDistrib)[1])->Arima)[2])->MA:=ma; + +Set nonLinFlt.new = { [[ + NameBlock addInp1 = LogAditiveEffect( + Real enabled = True, + Text segmentName = "Mixto_Submodel", + Serie observation, + Date First(observation), + Date Last(observation), + Set aditiveInputs = { SetOfNameBlock( + NameBlock addInp_X1 = [[ + //##Serie input + Serie _.serie = X1; + //##Si es falso no se usará este input + Real _.enabled = True; + //##Número de sigmas para la distribución a priori de los omitidos + Real _.mis_pri_sig_fac = 5; + //##Valor mínimo a priori para los omitidos + Real _.mis_pri_min = -1/0; + //##Valor máximo a priori para los omitidos + Real _.mis_pri_max = +1/0; + //##Valor inicial del parámetro + Real _.initValue = beta1; + //##Valor mínimo a priori del parámetro + Real _.minimum = -100; + //##Valor máximo a priori del parámetro + Real _.maximum = +100 + ]]) } + ) +]] }; + +Set nonLinFlt.old = + ((((bsr.parsed::SubModule)[1])->NoiseDistrib)[1])->NonLinFilters; + +Set ((((bsr.parsed::SubModule)[1])->NoiseDistrib)[1])->NonLinFilters := + nonLinFlt.new; + +BysMcmc::@Cycler cycler = BysMcmc::Bsr::Gibbs::GetCycler( + bsr.parsed, + BysMcmc::@McmcPartialHandler::Null(0), + bsr.config, + BysMcmc::@Notifier::Null(0)); + +BysMcmc::@Estimation estim = BysMcmc::@Estimation::New(cycler); + +Real estim::Run(0); + +/* */ Index: tol_tests/tol/Bugzilla/bug_001608/runMms.tol =================================================================== --- tol_tests/tol/Bugzilla/bug_001608/runMms.tol (revision b27f50548784bca303be5018af783b0f3f8a3f90) +++ tol_tests/tol/Bugzilla/bug_001608/runMms.tol (revision b27f50548784bca303be5018af783b0f3f8a3f90) @@ -0,0 +1,162 @@ +////////////////////////////////////////////////////////////////////////////// +// File: ejemplo1.tol +// Purpuse: ejemplo inventado de un modelo mixto +////////////////////////////////////////////////////////////////////////////// + +////////////////////////////////////////////////////////////////////////////// +#Require MMS; +////////////////////////////////////////////////////////////////////////////// + + + +////////////////////////////////////////////////////////////////////////////// +// 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 = ? + ]] +]]); + +// 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(?); +////////////////////////////////////////////////////////////////////////////// + +Real estimationBSR2::GetStrategy(?)::SaveIsolationFile("estimationBSR2.zip"); +