From bba434298be94e18314fef3fa032b84301c384f6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tomasz=20Wo=C5=BAniak?= Date: Thu, 6 Jun 2024 22:23:08 +1000 Subject: [PATCH 1/8] add small descriprtion #3 --- inst/varia/ilo_cubic_panel.R | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/inst/varia/ilo_cubic_panel.R b/inst/varia/ilo_cubic_panel.R index d2e7a39..edaa5b4 100644 --- a/inst/varia/ilo_cubic_panel.R +++ b/inst/varia/ilo_cubic_panel.R @@ -1,6 +1,5 @@ -# A first look at the data - +# Create ILO dataset from provided files library(dplyr) # this file contains a cubic panel dynamic dataset From 0249122ecb841a8cde9035832d8c0f1e6bdf588b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tomasz=20Wo=C5=BAniak?= Date: Thu, 6 Jun 2024 22:23:35 +1000 Subject: [PATCH 2/8] create ilo_conditional_forecast file #3 --- inst/varia/ilo_conditional_forecast.R | 26 ++++++++++++++++++++++++++ 1 file changed, 26 insertions(+) create mode 100644 inst/varia/ilo_conditional_forecast.R diff --git a/inst/varia/ilo_conditional_forecast.R b/inst/varia/ilo_conditional_forecast.R new file mode 100644 index 0000000..60106fd --- /dev/null +++ b/inst/varia/ilo_conditional_forecast.R @@ -0,0 +1,26 @@ + +# Create ILO conditional forecasts from provided files +library(dplyr) + +# this file contains a cubic panel dynamic dataset +all_cv <- read.csv("inst/varia/ilo_cubic_panel.csv") +colnames(all_cv) = c("year", "iso3code", "country", "UR", "EPR", "LFPR", "dgdp") + +# all variables all countries +data_cv <- all_cv %>% + filter(year >= 2024) %>% + select(year, iso3code, UR, EPR, LFPR, dgdp) + +# Create a list with the country data +countries = unique(data_cv$iso3code) +countries = countries[order(countries)] +ilo_conditional_forecast = list() +for (i in 1:length(countries)) { + ilo_conditional_forecast[[i]] <- data_cv %>% + filter(iso3code == countries[i]) %>% + select(UR, EPR, LFPR, dgdp) %>% + ts(start = 2024, frequency = 1) + names(ilo_conditional_forecast)[i] <- countries[i] +} + +save(ilo_conditional_forecast, file = "data/ilo_conditional_forecast.rda") From 4295fce8d8b60f3fc314b232eb80197b9f17633e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tomasz=20Wo=C5=BAniak?= Date: Thu, 6 Jun 2024 22:40:29 +1000 Subject: [PATCH 3/8] ilo_conditional_forecast data included #3 + data file + documentation --- R/ilo_conditional_forecast.R | 29 ++++++++++++++++++++++++ data/ilo_conditional_forecast.rda | Bin 0 -> 5956 bytes man/ilo_conditional_forecast.Rd | 36 ++++++++++++++++++++++++++++++ 3 files changed, 65 insertions(+) create mode 100644 R/ilo_conditional_forecast.R create mode 100644 data/ilo_conditional_forecast.rda create mode 100644 man/ilo_conditional_forecast.Rd diff --git a/R/ilo_conditional_forecast.R b/R/ilo_conditional_forecast.R new file mode 100644 index 0000000..a416393 --- /dev/null +++ b/R/ilo_conditional_forecast.R @@ -0,0 +1,29 @@ + +#' @title Data containing conditional projections for growth rate of GDP (dgdp) +#' for 189 United Nations countries from 2024 to 2029 +#' +#' @description For each of the countries a time series of 6 observations on +#' GDP growth rates (sgdp) #' formatted so they is provided to generate +#' conditional forecasts of labour market outcomes given the provided projected +#' paths of output. Last data update was implemented on 2024-05-11. +#' +#' @usage data(ilo_conditional_forecast) +#' +#' @format A list of 189 \code{ts} objects with time series of 6 observations +#' on 4 variables: +#' \describe{ +#' \item{UR}{unemployment rate - contains missing values} +#' \item{EPR}{annual employment rate - contains missing values} +#' \item{LFPR}{annual labour force participation rate - contains missing values} +#' \item{dgdp}{annual growth rate of gross domestic product - contains projected +#' values} +#' } +#' +#' @source +#' International Labour Organization. (2020). ILO modelled estimates database, +#' ILOSTAT [database]. Available from \url{https://ilostat.ilo.org/data/}. +#' +#' @examples +#' data(ilo_conditional_forecast) # upload the data +#' +"ilo_conditional_forecast" \ No newline at end of file diff --git a/data/ilo_conditional_forecast.rda b/data/ilo_conditional_forecast.rda new file mode 100644 index 0000000000000000000000000000000000000000..2449300d55b470514a98e644a787a9cfc8febaa9 GIT binary patch literal 5956 zcmXX}cRbYpACE|`3vt9ZJLAsC%!or|Zzr;IcJ|&QSsi7DvW1JRQ)brTvd7^fduL{y z5#e{&H~=bJ75=Cyy{^_`sHfTy^3?Zcs~a;Y;8Ua1kB29HZ^ z(DjZd0ChD*)~K1RYo*jQAOYc9QiI?|P7XTsXQ zHF&Yx|4`6M({BKUe3m|u6TDc{Y$VJ(Hs&xO=dihLGFCF2B^Uf7jpR!J+=stzh~PIqLZzK2^Z4 z*q<8!a!fk5w7-%mz!1CYZttvn%lNxPwJNXAlhxUsb5K}Q*kSA`=} zawTV3V>fEAdUjYq?$Gh(Nr>q38Dd^xC8KjTc`7{iZ+qOwy1^?kIY$3^JY4qx8DZFN`v;XT5 zPZhhaB+Zr4d=oVp?aj*il@;#IR$+?TmuEjV!q3D^_$mH9_cLxxD_~68I=C_8F7K)q zwDiTT5k(1~H7jG&vZRXcmXW79)jj`Z=P_NLdww?tUWhE}q}RN=7D{apfTX#aP9ZEV zboTn=%RAwEZA(kF^HqFM+?Exav;@+!y*6;~Gn+2!ZMR_1^nP#*oV3s0Y`7`}ww<%S zPK7#ma(Kiu`WR*ityI%; z+S6*7R}3vyF({W(kQqUH4?nZhb|HML^e~4_?MmsOCc7-w#|2>h;+h}8LiXT=qx96LTqT%g^cgF^lj zw;1u;&KtX-nHy(SS25y>HejCoB1gZc3Gqs5E*NV8Px+QG>~D5sg^s!=0gE0&%ORZ6 zmKu-jBGT}e*Fs1QL_u_VA#BsP-UaQ^0k1H1BHO$RACC~bUzB4jn69J)*nGZA68fRh zWf?OwMuncalHMV0-Is{A-&+hK4{^?xWEVm@dR=Ex5+T97*`h{K$MT~x>@9-#hR!t! z>kxN$X-Ma_h^NT%BFMuVHxyB060WoLy4!WCUxapPzWtY7+`PW-p>+z6b$!Bm!b;eP zJ0V78@6pG}k*ZW6CL-+-nfXOgmc7jdmy!q;X8w^`o0{RzO2zrMO7VKX;BtJJG_!!@ zl=0c#2@pe;7P83{coB4CHZZ(bF994>6Bw8{B@sS$X|bT%1|W01!_)wr?F2nX6WNzs zDz`5e#GaZsZD<+W1MjA%j7BbPt;lxk{Uhfj-K_o?d)V8L=Ik1B=Zf6{GXO07Z)kxl zx~eI@ATu+x_egD+llO-j{-?_F!B)GMz%gc_uJ&)FPU=%%`)1}ZNEKif#Z znjM=+Giw{R@qICYz=;DYJ|7eXed9VB%z(Xroi7{vS=$HRMYtk~og|+JR$ZTeRQTu1 zlJm%TeD|5Go3>fe~y1S99ae?34DhN_eqInpIVVRWI`n8Dk9|6pxk-onCQ$# z^60NiiEA=2n;1C!xONP?>=5ybDLCjz;M=RV>$DvWUa^iQZZ~F*0A0dOtZ>_N#0AYq zQvk!1da@l56Q{%BV3fVkbyww-E~#-M`IGx5hYXToC?OqYn8CD{HzfsS^rloRYC5)R5aRZ#De+??lyq z+@4r?gh^}-*nPv?`@O&~^K zdO2G3W>x$OmFTRys!~xD;V00lxB(m`4msl#RbW`mImsZ94RlwL#$v^DyJOrH`}#}8 zC2`nO6z65Q)H1koR(pP$vhI-m>YX6cl9Vrwj(w{-<9SazWxNZVfrnZIM284_w|ra9 zAVlRSI&>w_$0D&Lc(_wXM>0UfF|iqnruZV3&K#)RkZVnV zcID3lGrZBy7vScgVj*wC{{UC67b$O)1mosL0tZs8S7|PewC~GN*DK;r5=|>9D&d4R zr;5RwbI&Azr13W1t-w&=B|5c)diFr}RMMs4`3iMET7JsRdsyS1VwHAxG0$aD6fjOk#Oqwz zcG;kt{j?2f8MGZP);IN0uh=&kh&mn;xkitThWvzAUGRw-`4p5|`V{Eg`G(aV(9^<$ zQA#rp1?$=B*2MTK8mN5fIkU8d>MYSE z{lkg36INR70DlFw#NlHJa-jIyK)sX`z|B_@vC9eAeo)x{nsZ23go69xj?U|NNTS^P z>WtsBNkxCpxQ6xnl-UY9-vA*Z^t3e%V zLf?ZB%7OTn>1@)jLJHpmfb^91OZsmOKU$>0cbb)+O7Z3~#FH1n0NYnW)12H&A|fFz z&7nsk^_Tk$C}_*d-+{M&`y2G$VlOM&Xspr%FD{2GD`~6*q|RU87j^E+(50Fh36ym3 zV%mUI)1^kM$i=Ls)7=stUw^XA=P|&FEPn(Hokp#S5hVC8h-u?vgV3&ud0v2B7_{WY ztBl(zH+%t4H;0ggBjjJ>1iDBcy|@Ps8YB!svh7#iW`=ZDn$3(y4iFKCz?cL-KJV@}$K&b=ioo)GAiWm;v32c4|0wo`gTxqW099^OoP zu%FC3`sR@ez5s22w;m)ob0N}abV)>6|Cs!1mPOPI`i`lrij?9ga~i85zA%xK+oEQK zH;}o73AsWRkWC&QMVQE%3<9*e){D%IFTFDX8&5o9IpP+6L=va6BD0jc4z5eHrJmu^ z;kCI0!+(t131B9ujj|DJrywM#@nY;TL3~8cQ|bigE;mLx^`NU_eOd0^qU`W{d;$9$ z%Cq4Lfp~WtI3LlRNYhd*p4d@Z;20A!u4F$L&txu=v;HH|_LpI|FX z|D}A-ZxHFBnn_%!UQMBm<*#AS+scwj@8vE5b7-+g z4qYA~+Ms@XxrtjNERFHK)e?-wS%g^b9IKW8=;q6->!|x|4J&_^gV^_;4R1jL1ov|i%WYm=m0r}u(ms|zJZVFrd=kF+~8`#vK-q+I{>BlJw7djmt@i*&Zu zp`F($r?Pu8#b`9Esz+Y$9YVDM+XJz|Q)ory2glN>mHbj&>q%5r{BF6e8wGxG=8b=H z1cgpHLNHw{_vcM0lhK;>u7Rg{?5(v`A-#)%G&0Wnl)BAQSp|$65JRU9&Dv+OhE9=v zbdk#M4jfxL2*-M2mqwIEwF4m!L^r;H)rK)Lh5x0(GoX!A@W&1jEe?MeQ_*J6d zK9@Mhp5|uKP8YhjP%t1k%sWKOMr59c2*K@vR-N;zuU3*7(*Ys~IAxFE2!x(90`rrR zU`US*_Bdq9ocHqBwYE})Ti8k(s)gN#h!G)K?u-m~q>LqQncABhYFK3~e+y{ARUe2_ zzK1l)j(b6c&}DeYUP`O0-q$tW^MB7S42ae6=ejvWj z_#CYp16M#mixvb0kW>t~i<<4YLQP9D=QuhuG9!pgW5Q8d$bq<*zC;Sxt*t2O3=>_YryAEqv91c~qHXXP0MC$yP9`@n{mK_&mLeE-p+707P!y-= zVmG*|H7?siB}V7W~uxew8SP8qI3JK`(%bC#C{$p;_rXQP&k zUHNrrgOP%#=&eA0W#jEydFJo{t!>DKpqwhHk+JXW&{(n4U1vAu)>$> zn{IWKkIUmINH!tp{)U($)MET^y}_EILH<8DqB$)U$3B22VA}+l5O=^uN+AvN_Zxbn zbW5o(2>w9AeIWx(Q>#+xv04qmJ^N<*tPmu#)`&xZUX5Oa zY%a}eS8dke4{EHLS9!=)?-^Aj z0tuQAS$}ufvI@`JP24#dx_9Hpedygms!FnxTK@;%HdJ>bf_JvNUzDR3_GeGGEpF`+ zF~=wM&K0)2+IXuhv?hnJbh@|vGebIu^E%I&^2C>oOPW64Zu6$@Fa9&|)m*bH#^OOTI52Cq}~gsw~%C(7rVv0Gi&KWc+%BYOXmIj1F?1Tod$y?(-)&Y@X5!}LPWK}(MWb$L#K6e_ zc2{S;10A&{R4rujMte2w4!JX=IXu^h*_<@`2d`+!DZ4Nu;xWXsP)+UbPTo>s?viN& zPw=fngZ2(uLm=jfZZ9%NPw7>Tv2+(Cg@8bk)^ec5$SkcFcl1Izj)JxiMG zINHU3LWCfw5PRv_vBS}74uSx8!gzudMg={0Gu_JN8VY$1@<(zMF%3xgX8Pn4}<=nwPD$51;& zN=1sxzH8oz!YM^^@+?n^8OG^@CWMGj7N(O`q;AC)%)6Dn|h7svVsuzV0ByO;wnr!`~l{Kb?*#D=ppc|$aEt3$MEUl=9A2G6i%Kh2A?&L-Dx z;9Ppw?TsEzk9VsW|0uJ6+TtcrwwCFIE3X{aK2lJfU>-4?muRX}dAM~%a^zYjtrKqc z3^jDR605sG&|n?qpr;URZ|bb^IelkMb;T`dHk@# zan!-lNVld?c~%)y#Dg2xZK)|!e%M`>_*CLYX~`v`Pr_H}DgCC+Q(UScEZp;HiHc;F z4=lN=uk8Bw@v>jSDpsA0gN04zv&u(~M?Phr$JyYuVdJQ5V-Jk8fqhHwL zT#f1p6`vTvkk`h=`2-#gZkA5i;)*A}esca)F|K`qnpz{2NFn~2|MTGbbN4uh4h+H0JOrEB_`gWo zMoKpC+%^y2!gP@ZCgEC%=C|qafxM;u_a4T72Vg2t8lh%3wN3HgO?|W89Bl^Odl>PZ zCp1O^$2^gLP5j-&MoE})xLGQz^#^>*!wvn9>i%OT_MxJo-P(fy7nln4Q;GVEw0l;S zZ(UiD?GjrQxN?eat;R>}A7zwOuC(pTaHmOF05{C4I1AvflBmK2Rd5j9RJX4=vo*!p zn{N5Q!uQE8=;=Uz2Dpw;F42$NjhtJ`SzT&uaDiWZH~VUiHl=O5OWphJkC^bMNw~Z3 zbMohL_g&*Jr=3b+-F1$Nf=r Date: Sat, 8 Jun 2024 23:03:33 +1000 Subject: [PATCH 4/8] develop cpp function for RNG from conditional mvnorm #8 --- inst/include/bvarPANELs_RcppExports.h | 21 +++++++++++++++ src/RcppExports.cpp | 39 +++++++++++++++++++++++++++ src/forecast_panel.cpp | 34 +++++++++++++++++++++++ src/forecast_panel.h | 7 +++++ 4 files changed, 101 insertions(+) diff --git a/inst/include/bvarPANELs_RcppExports.h b/inst/include/bvarPANELs_RcppExports.h index f5bfa69..86e0ff0 100644 --- a/inst/include/bvarPANELs_RcppExports.h +++ b/inst/include/bvarPANELs_RcppExports.h @@ -46,6 +46,27 @@ namespace bvarPANELs { return Rcpp::as(rcpp_result_gen); } + inline arma::vec mvnrnd_cond(arma::vec x, arma::vec mu, arma::mat Sigma) { + typedef SEXP(*Ptr_mvnrnd_cond)(SEXP,SEXP,SEXP); + static Ptr_mvnrnd_cond p_mvnrnd_cond = NULL; + if (p_mvnrnd_cond == NULL) { + validateSignature("arma::vec(*mvnrnd_cond)(arma::vec,arma::vec,arma::mat)"); + p_mvnrnd_cond = (Ptr_mvnrnd_cond)R_GetCCallable("bvarPANELs", "_bvarPANELs_mvnrnd_cond"); + } + RObject rcpp_result_gen; + { + RNGScope RCPP_rngScope_gen; + rcpp_result_gen = p_mvnrnd_cond(Shield(Rcpp::wrap(x)), Shield(Rcpp::wrap(mu)), Shield(Rcpp::wrap(Sigma))); + } + if (rcpp_result_gen.inherits("interrupted-error")) + throw Rcpp::internal::InterruptedException(); + if (Rcpp::internal::isLongjumpSentinel(rcpp_result_gen)) + throw Rcpp::LongjumpException(rcpp_result_gen); + if (rcpp_result_gen.inherits("try-error")) + throw Rcpp::exception(Rcpp::as(rcpp_result_gen).c_str()); + return Rcpp::as(rcpp_result_gen); + } + } #endif // RCPP_bvarPANELs_RCPPEXPORTS_H_GEN_ diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index 47765c9..d48e1cc 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -69,6 +69,42 @@ RcppExport SEXP _bvarPANELs_forecast_bvarPANEL(SEXP posterior_A_c_cppSEXP, SEXP UNPROTECT(1); return rcpp_result_gen; } +// mvnrnd_cond +arma::vec mvnrnd_cond(arma::vec x, arma::vec mu, arma::mat Sigma); +static SEXP _bvarPANELs_mvnrnd_cond_try(SEXP xSEXP, SEXP muSEXP, SEXP SigmaSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::traits::input_parameter< arma::vec >::type x(xSEXP); + Rcpp::traits::input_parameter< arma::vec >::type mu(muSEXP); + Rcpp::traits::input_parameter< arma::mat >::type Sigma(SigmaSEXP); + rcpp_result_gen = Rcpp::wrap(mvnrnd_cond(x, mu, Sigma)); + return rcpp_result_gen; +END_RCPP_RETURN_ERROR +} +RcppExport SEXP _bvarPANELs_mvnrnd_cond(SEXP xSEXP, SEXP muSEXP, SEXP SigmaSEXP) { + SEXP rcpp_result_gen; + { + Rcpp::RNGScope rcpp_rngScope_gen; + rcpp_result_gen = PROTECT(_bvarPANELs_mvnrnd_cond_try(xSEXP, muSEXP, SigmaSEXP)); + } + Rboolean rcpp_isInterrupt_gen = Rf_inherits(rcpp_result_gen, "interrupted-error"); + if (rcpp_isInterrupt_gen) { + UNPROTECT(1); + Rf_onintr(); + } + bool rcpp_isLongjump_gen = Rcpp::internal::isLongjumpSentinel(rcpp_result_gen); + if (rcpp_isLongjump_gen) { + Rcpp::internal::resumeJump(rcpp_result_gen); + } + Rboolean rcpp_isError_gen = Rf_inherits(rcpp_result_gen, "try-error"); + if (rcpp_isError_gen) { + SEXP rcpp_msgSEXP_gen = Rf_asChar(rcpp_result_gen); + UNPROTECT(1); + Rf_error("%s", CHAR(rcpp_msgSEXP_gen)); + } + UNPROTECT(1); + return rcpp_result_gen; +} // rmniw1 arma::field rmniw1(const arma::mat& A, const arma::mat& V, const arma::mat& S, const double& nu); RcppExport SEXP _bvarPANELs_rmniw1(SEXP ASEXP, SEXP VSEXP, SEXP SSEXP, SEXP nuSEXP) { @@ -225,6 +261,7 @@ static int _bvarPANELs_RcppExport_validate(const char* sig) { static std::set signatures; if (signatures.empty()) { signatures.insert("Rcpp::List(*forecast_bvarPANEL)(arma::field&,arma::field&,Rcpp::List&,const int)"); + signatures.insert("arma::vec(*mvnrnd_cond)(arma::vec,arma::vec,arma::mat)"); } return signatures.find(sig) != signatures.end(); } @@ -232,6 +269,7 @@ static int _bvarPANELs_RcppExport_validate(const char* sig) { // registerCCallable (register entry points for exported C++ functions) RcppExport SEXP _bvarPANELs_RcppExport_registerCCallable() { R_RegisterCCallable("bvarPANELs", "_bvarPANELs_forecast_bvarPANEL", (DL_FUNC)_bvarPANELs_forecast_bvarPANEL_try); + R_RegisterCCallable("bvarPANELs", "_bvarPANELs_mvnrnd_cond", (DL_FUNC)_bvarPANELs_mvnrnd_cond_try); R_RegisterCCallable("bvarPANELs", "_bvarPANELs_RcppExport_validate", (DL_FUNC)_bvarPANELs_RcppExport_validate); return R_NilValue; } @@ -239,6 +277,7 @@ RcppExport SEXP _bvarPANELs_RcppExport_registerCCallable() { static const R_CallMethodDef CallEntries[] = { {"_bvarPANELs_bvarPANEL", (DL_FUNC) &_bvarPANELs_bvarPANEL, 8}, {"_bvarPANELs_forecast_bvarPANEL", (DL_FUNC) &_bvarPANELs_forecast_bvarPANEL, 4}, + {"_bvarPANELs_mvnrnd_cond", (DL_FUNC) &_bvarPANELs_mvnrnd_cond, 3}, {"_bvarPANELs_rmniw1", (DL_FUNC) &_bvarPANELs_rmniw1, 4}, {"_bvarPANELs_sample_m", (DL_FUNC) &_bvarPANELs_sample_m, 5}, {"_bvarPANELs_sample_w", (DL_FUNC) &_bvarPANELs_sample_w, 2}, diff --git a/src/forecast_panel.cpp b/src/forecast_panel.cpp index 346e7a4..be6e2b5 100644 --- a/src/forecast_panel.cpp +++ b/src/forecast_panel.cpp @@ -56,3 +56,37 @@ Rcpp::List forecast_bvarPANEL ( _["forecasts_cpp"] = forecasts ); } // END forecast_bsvarPANEL + + + +// [[Rcpp::interfaces(cpp)]] +// [[Rcpp::export]] +arma::vec mvnrnd_cond ( + arma::vec x, // Nx1 with NAs or without + arma::vec mu, // Nx1 mean vector + arma::mat Sigma // NxN covariance matrix +) { + int N = x.n_elem; + uvec ind = find_finite(x); + uvec ind_nan = find_nan(x); + mat aj = eye(N, N); + + vec x2 = x(ind); + + vec mu1 = mu(ind_nan); + vec mu2 = mu(ind); + mat Sigma11 = Sigma(ind_nan, ind_nan); + mat Sigma12 = Sigma(ind_nan, ind); + mat Sigma22 = Sigma(ind, ind); + mat Sigma22_inv = inv_sympd(Sigma22); + + vec mu_cond = mu1 + Sigma12 * Sigma22_inv * (x2 - mu2); + mat Sigma_cond = Sigma11 - Sigma12 * Sigma22_inv * Sigma12.t(); + + vec draw = mvnrnd( mu_cond, Sigma_cond); + + vec out = aj.cols(ind_nan) * draw + aj.cols(ind) * x2; + return out; +} // END mvnrnd_cond + + diff --git a/src/forecast_panel.h b/src/forecast_panel.h index 44d6290..95e21ca 100644 --- a/src/forecast_panel.h +++ b/src/forecast_panel.h @@ -15,4 +15,11 @@ Rcpp::List forecast_bvarPANEL ( ); +arma::vec mvnrnd_cond ( + arma::vec x, // Nx1 with NAs or without + arma::vec mu, // Nx1 mean vector + arma::mat Sigma // NxN covariance matrix +); + + #endif // _FORECAST_PANEL_H_ From b087056cd405d4b862a649673694e1c8dc68e0d1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tomasz=20Wo=C5=BAniak?= Date: Mon, 10 Jun 2024 18:34:13 +1000 Subject: [PATCH 5/8] developped forecast_conditional_bvarPANEL #8 --- inst/include/bvarPANELs_RcppExports.h | 21 ++++++++++ src/RcppExports.cpp | 41 +++++++++++++++++++ src/forecast_panel.cpp | 57 +++++++++++++++++++++++++++ src/forecast_panel.h | 9 +++++ 4 files changed, 128 insertions(+) diff --git a/inst/include/bvarPANELs_RcppExports.h b/inst/include/bvarPANELs_RcppExports.h index 86e0ff0..7b78d45 100644 --- a/inst/include/bvarPANELs_RcppExports.h +++ b/inst/include/bvarPANELs_RcppExports.h @@ -67,6 +67,27 @@ namespace bvarPANELs { return Rcpp::as(rcpp_result_gen); } + inline Rcpp::List forecast_conditional_bvarPANEL(arma::field& posterior_A_c_cpp, arma::field& posterior_Sigma_c_cpp, Rcpp::List& X_c, Rcpp::List& cond_forecasts, const int horizon) { + typedef SEXP(*Ptr_forecast_conditional_bvarPANEL)(SEXP,SEXP,SEXP,SEXP,SEXP); + static Ptr_forecast_conditional_bvarPANEL p_forecast_conditional_bvarPANEL = NULL; + if (p_forecast_conditional_bvarPANEL == NULL) { + validateSignature("Rcpp::List(*forecast_conditional_bvarPANEL)(arma::field&,arma::field&,Rcpp::List&,Rcpp::List&,const int)"); + p_forecast_conditional_bvarPANEL = (Ptr_forecast_conditional_bvarPANEL)R_GetCCallable("bvarPANELs", "_bvarPANELs_forecast_conditional_bvarPANEL"); + } + RObject rcpp_result_gen; + { + RNGScope RCPP_rngScope_gen; + rcpp_result_gen = p_forecast_conditional_bvarPANEL(Shield(Rcpp::wrap(posterior_A_c_cpp)), Shield(Rcpp::wrap(posterior_Sigma_c_cpp)), Shield(Rcpp::wrap(X_c)), Shield(Rcpp::wrap(cond_forecasts)), Shield(Rcpp::wrap(horizon))); + } + if (rcpp_result_gen.inherits("interrupted-error")) + throw Rcpp::internal::InterruptedException(); + if (Rcpp::internal::isLongjumpSentinel(rcpp_result_gen)) + throw Rcpp::LongjumpException(rcpp_result_gen); + if (rcpp_result_gen.inherits("try-error")) + throw Rcpp::exception(Rcpp::as(rcpp_result_gen).c_str()); + return Rcpp::as(rcpp_result_gen); + } + } #endif // RCPP_bvarPANELs_RCPPEXPORTS_H_GEN_ diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index d48e1cc..a8865ca 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -105,6 +105,44 @@ RcppExport SEXP _bvarPANELs_mvnrnd_cond(SEXP xSEXP, SEXP muSEXP, SEXP SigmaSEXP) UNPROTECT(1); return rcpp_result_gen; } +// forecast_conditional_bvarPANEL +Rcpp::List forecast_conditional_bvarPANEL(arma::field& posterior_A_c_cpp, arma::field& posterior_Sigma_c_cpp, Rcpp::List& X_c, Rcpp::List& cond_forecasts, const int horizon); +static SEXP _bvarPANELs_forecast_conditional_bvarPANEL_try(SEXP posterior_A_c_cppSEXP, SEXP posterior_Sigma_c_cppSEXP, SEXP X_cSEXP, SEXP cond_forecastsSEXP, SEXP horizonSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::traits::input_parameter< arma::field& >::type posterior_A_c_cpp(posterior_A_c_cppSEXP); + Rcpp::traits::input_parameter< arma::field& >::type posterior_Sigma_c_cpp(posterior_Sigma_c_cppSEXP); + Rcpp::traits::input_parameter< Rcpp::List& >::type X_c(X_cSEXP); + Rcpp::traits::input_parameter< Rcpp::List& >::type cond_forecasts(cond_forecastsSEXP); + Rcpp::traits::input_parameter< const int >::type horizon(horizonSEXP); + rcpp_result_gen = Rcpp::wrap(forecast_conditional_bvarPANEL(posterior_A_c_cpp, posterior_Sigma_c_cpp, X_c, cond_forecasts, horizon)); + return rcpp_result_gen; +END_RCPP_RETURN_ERROR +} +RcppExport SEXP _bvarPANELs_forecast_conditional_bvarPANEL(SEXP posterior_A_c_cppSEXP, SEXP posterior_Sigma_c_cppSEXP, SEXP X_cSEXP, SEXP cond_forecastsSEXP, SEXP horizonSEXP) { + SEXP rcpp_result_gen; + { + Rcpp::RNGScope rcpp_rngScope_gen; + rcpp_result_gen = PROTECT(_bvarPANELs_forecast_conditional_bvarPANEL_try(posterior_A_c_cppSEXP, posterior_Sigma_c_cppSEXP, X_cSEXP, cond_forecastsSEXP, horizonSEXP)); + } + Rboolean rcpp_isInterrupt_gen = Rf_inherits(rcpp_result_gen, "interrupted-error"); + if (rcpp_isInterrupt_gen) { + UNPROTECT(1); + Rf_onintr(); + } + bool rcpp_isLongjump_gen = Rcpp::internal::isLongjumpSentinel(rcpp_result_gen); + if (rcpp_isLongjump_gen) { + Rcpp::internal::resumeJump(rcpp_result_gen); + } + Rboolean rcpp_isError_gen = Rf_inherits(rcpp_result_gen, "try-error"); + if (rcpp_isError_gen) { + SEXP rcpp_msgSEXP_gen = Rf_asChar(rcpp_result_gen); + UNPROTECT(1); + Rf_error("%s", CHAR(rcpp_msgSEXP_gen)); + } + UNPROTECT(1); + return rcpp_result_gen; +} // rmniw1 arma::field rmniw1(const arma::mat& A, const arma::mat& V, const arma::mat& S, const double& nu); RcppExport SEXP _bvarPANELs_rmniw1(SEXP ASEXP, SEXP VSEXP, SEXP SSEXP, SEXP nuSEXP) { @@ -262,6 +300,7 @@ static int _bvarPANELs_RcppExport_validate(const char* sig) { if (signatures.empty()) { signatures.insert("Rcpp::List(*forecast_bvarPANEL)(arma::field&,arma::field&,Rcpp::List&,const int)"); signatures.insert("arma::vec(*mvnrnd_cond)(arma::vec,arma::vec,arma::mat)"); + signatures.insert("Rcpp::List(*forecast_conditional_bvarPANEL)(arma::field&,arma::field&,Rcpp::List&,Rcpp::List&,const int)"); } return signatures.find(sig) != signatures.end(); } @@ -270,6 +309,7 @@ static int _bvarPANELs_RcppExport_validate(const char* sig) { RcppExport SEXP _bvarPANELs_RcppExport_registerCCallable() { R_RegisterCCallable("bvarPANELs", "_bvarPANELs_forecast_bvarPANEL", (DL_FUNC)_bvarPANELs_forecast_bvarPANEL_try); R_RegisterCCallable("bvarPANELs", "_bvarPANELs_mvnrnd_cond", (DL_FUNC)_bvarPANELs_mvnrnd_cond_try); + R_RegisterCCallable("bvarPANELs", "_bvarPANELs_forecast_conditional_bvarPANEL", (DL_FUNC)_bvarPANELs_forecast_conditional_bvarPANEL_try); R_RegisterCCallable("bvarPANELs", "_bvarPANELs_RcppExport_validate", (DL_FUNC)_bvarPANELs_RcppExport_validate); return R_NilValue; } @@ -278,6 +318,7 @@ static const R_CallMethodDef CallEntries[] = { {"_bvarPANELs_bvarPANEL", (DL_FUNC) &_bvarPANELs_bvarPANEL, 8}, {"_bvarPANELs_forecast_bvarPANEL", (DL_FUNC) &_bvarPANELs_forecast_bvarPANEL, 4}, {"_bvarPANELs_mvnrnd_cond", (DL_FUNC) &_bvarPANELs_mvnrnd_cond, 3}, + {"_bvarPANELs_forecast_conditional_bvarPANEL", (DL_FUNC) &_bvarPANELs_forecast_conditional_bvarPANEL, 5}, {"_bvarPANELs_rmniw1", (DL_FUNC) &_bvarPANELs_rmniw1, 4}, {"_bvarPANELs_sample_m", (DL_FUNC) &_bvarPANELs_sample_m, 5}, {"_bvarPANELs_sample_w", (DL_FUNC) &_bvarPANELs_sample_w, 2}, diff --git a/src/forecast_panel.cpp b/src/forecast_panel.cpp index be6e2b5..e1c87e5 100644 --- a/src/forecast_panel.cpp +++ b/src/forecast_panel.cpp @@ -90,3 +90,60 @@ arma::vec mvnrnd_cond ( } // END mvnrnd_cond + + +// [[Rcpp::interfaces(cpp)]] +// [[Rcpp::export]] +Rcpp::List forecast_conditional_bvarPANEL ( + arma::field& posterior_A_c_cpp, // (S)(K, N, C) + arma::field& posterior_Sigma_c_cpp, // (S)(N, N, C) + Rcpp::List& X_c, // (C)(T_c, K) + Rcpp::List& cond_forecasts, // (C)(horizon, N) + const int horizon +) { + + const int N = posterior_A_c_cpp(0).n_cols; + const int S = posterior_A_c_cpp.n_elem; + const int K = posterior_A_c_cpp(0).n_rows; + const int C = posterior_A_c_cpp(0).n_slices; + const int p = (K - 1) / N; + + field forecasts(C); // of (horizon, N, S) cubes + rowvec one(1, fill::value(1)); + + for (int c=0; c(X_c[c]); + mat cond_fc = as(cond_forecasts[c]); + rowvec x_t = XXcc.tail_rows(1).cols(0, K - 2); + cube forecasts_c(horizon, N, S); + + for (int s=0; s& posterior_A_c_cpp, // (S)(K, N, C) + arma::field& posterior_Sigma_c_cpp, // (S)(N, N, C) + Rcpp::List& X_c, // (C)(T_c, K) + Rcpp::List& cond_forecasts, // (C)(horizon, N) + const int horizon +); + + #endif // _FORECAST_PANEL_H_ From 74b0dcb60507b6c4f40dc60c0f9c370fe4aea6a0 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tomasz=20Wo=C5=BAniak?= Date: Mon, 10 Jun 2024 19:26:55 +1000 Subject: [PATCH 6/8] expand arguments in generic forecast #8 --- R/forecast.R | 14 ++++++++++++-- man/forecast.PosteriorBVARPANEL.Rd | 14 +++++++++++++- 2 files changed, 25 insertions(+), 3 deletions(-) diff --git a/R/forecast.R b/R/forecast.R index fb5eb6a..9e47d64 100644 --- a/R/forecast.R +++ b/R/forecast.R @@ -11,6 +11,12 @@ #' @param horizon a positive integer, specifying the forecasting horizon. #' @param exogenous_forecast not used here ATM; included for compatibility with #' generic \code{forecast}. +#' @param conditional_forecast a list of length \code{C} containing +#' \code{horizon x N} matrices or \code{ts} objects with of the same dimensions +#' with forecasted values for selected variables. These objects should only +#' contain \code{numeric} or \code{NA} values. The entries with \code{NA} values +#' correspond to the values that are forecasted conditionally on the realisations +#' provided as \code{numeric} values. #' #' @return A list of class \code{PanelForecasts} containing the #' draws from the predictive density and data. The output list includes element: @@ -41,8 +47,12 @@ #' forecast(horizon = 2) -> predictive #' #' @export -forecast.PosteriorBVARPANEL = function(posterior, horizon = 1, exogenous_forecast = NULL) { - +forecast.PosteriorBVARPANEL = function( + posterior, + horizon = 1, + exogenous_forecast = NULL, + conditional_forecast = NULL +) { posterior_A_c_cpp = posterior$posterior$A_c_cpp posterior_Sigma_c_cpp = posterior$posterior$Sigma_c_cpp diff --git a/man/forecast.PosteriorBVARPANEL.Rd b/man/forecast.PosteriorBVARPANEL.Rd index b286fc6..cdf9f80 100644 --- a/man/forecast.PosteriorBVARPANEL.Rd +++ b/man/forecast.PosteriorBVARPANEL.Rd @@ -4,7 +4,12 @@ \alias{forecast.PosteriorBVARPANEL} \title{Forecasting using Hierarchical Pannel Vector Autoregressions} \usage{ -\method{forecast}{PosteriorBVARPANEL}(posterior, horizon = 1, exogenous_forecast = NULL) +\method{forecast}{PosteriorBVARPANEL}( + posterior, + horizon = 1, + exogenous_forecast = NULL, + conditional_forecast = NULL +) } \arguments{ \item{posterior}{posterior estimation outcome - an object of class @@ -14,6 +19,13 @@ PosteriorBVARPANEL obtained by running the \code{estimate} function.} \item{exogenous_forecast}{not used here ATM; included for compatibility with generic \code{forecast}.} + +\item{conditional_forecast}{a list of length \code{C} containing +\code{horizon x N} matrices or \code{ts} objects with of the same dimensions +with forecasted values for selected variables. These objects should only +contain \code{numeric} or \code{NA} values. The entries with \code{NA} values +correspond to the values that are forecasted conditionally on the realisations +provided as \code{numeric} values.} } \value{ A list of class \code{PanelForecasts} containing the From 0563e9db4702b8b1c8c8ef586be39fc1b5771c48 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tomasz=20Wo=C5=BAniak?= Date: Mon, 10 Jun 2024 19:59:59 +1000 Subject: [PATCH 7/8] expanded forecast.PosteriorBVARPANEL method to include conditional forecasting #8 --- R/forecast.R | 62 ++++++++++++++++++++++++++---- man/forecast.PosteriorBVARPANEL.Rd | 28 +++++++++++--- 2 files changed, 78 insertions(+), 12 deletions(-) diff --git a/R/forecast.R b/R/forecast.R index 9e47d64..9063b6f 100644 --- a/R/forecast.R +++ b/R/forecast.R @@ -12,11 +12,10 @@ #' @param exogenous_forecast not used here ATM; included for compatibility with #' generic \code{forecast}. #' @param conditional_forecast a list of length \code{C} containing -#' \code{horizon x N} matrices or \code{ts} objects with of the same dimensions -#' with forecasted values for selected variables. These objects should only -#' contain \code{numeric} or \code{NA} values. The entries with \code{NA} values -#' correspond to the values that are forecasted conditionally on the realisations -#' provided as \code{numeric} values. +#' \code{horizon x N} matrices with forecasted values for selected variables. +#' These matrices should only contain \code{numeric} or \code{NA} values. The +#' entries with \code{NA} values correspond to the values that are forecasted +#' conditionally on the realisations provided as \code{numeric} values. #' #' @return A list of class \code{PanelForecasts} containing the #' draws from the predictive density and data. The output list includes element: @@ -46,6 +45,25 @@ #' estimate(S = 20) |> #' forecast(horizon = 2) -> predictive #' +#' # conditional forecasting 6 years ahead conditioning on +#' # provided future values for the Gross Domestic Product +#' # growth rate +#' ############################################################ +#' #' data(ilo_conditional_forecast) # load the conditional forecasts of dgdp +#' predictive = forecast(posterior, 6, conditional_forecast = ilo_conditional_forecast) +#' +#' # workflow with the pipe |> +#' ############################################################ +#' set.seed(123) +#' ilo_cubic_panel |> +#' specify_bvarPANEL$new() |> +#' estimate(S = 10) |> +#' estimate(S = 20) |> +#' forecast( +#' horizon = 6, +#' conditional_forecast = ilo_conditional_forecast +#' ) -> predictive +#' #' @export forecast.PosteriorBVARPANEL = function( posterior, @@ -58,10 +76,40 @@ forecast.PosteriorBVARPANEL = function( posterior_Sigma_c_cpp = posterior$posterior$Sigma_c_cpp X_c = posterior$last_draw$data_matrices$X Y_c = posterior$last_draw$data_matrices$Y + N = dim(Y_c[[1]])[2] - fore = .Call(`_bvarPANELs_forecast_bvarPANEL`, posterior_A_c_cpp, posterior_Sigma_c_cpp, X_c, horizon) + do_conditional_forecasting = !is.null(conditional_forecast) + + if (!do_conditional_forecasting) { + # perform forecasting + fore = .Call(`_bvarPANELs_forecast_bvarPANEL`, + posterior_A_c_cpp, posterior_Sigma_c_cpp, X_c, horizon) + } else { + stopifnot("Argument conditional_forecast must be a list with the same countries + as in the provided data." + = is.list(conditional_forecast) & length(conditional_forecast) == length(Y_c) + ) + stopifnot("Argument conditional_forecast must be a list with the same countries + as in the provided data." + = all(names(Y_c) == names(conditional_forecast)) + ) + stopifnot("Argument conditional_forecast must be a list with matrices with numeric values." + = all(sapply(conditional_forecast, function(x) is.matrix(x) & is.numeric(x))) + ) + stopifnot("All the matrices provided in argument conditional_forecast must have + the same number of rows equal to the value of argument horizon." + = unique(sapply(conditional_forecast, function(x) nrow(x) )) == horizon + ) + stopifnot("All the matrices provided in argument conditional_forecast must have + the same number of columns equal to the number of columns in the used data." + = unique(sapply(conditional_forecast, function(x) ncol(x) )) == N + ) + + # perform conditional forecasting + fore = .Call(`_bvarPANELs_forecast_conditional_bvarPANEL`, + posterior_A_c_cpp, posterior_Sigma_c_cpp, X_c, conditional_forecast, horizon) + } - N = dim(Y_c[[1]])[2] S = dim(posterior_A_c_cpp)[1] C = length(Y_c) forecasts = array(NA, c(horizon, N, S, C)) diff --git a/man/forecast.PosteriorBVARPANEL.Rd b/man/forecast.PosteriorBVARPANEL.Rd index cdf9f80..6c3fb43 100644 --- a/man/forecast.PosteriorBVARPANEL.Rd +++ b/man/forecast.PosteriorBVARPANEL.Rd @@ -21,11 +21,10 @@ PosteriorBVARPANEL obtained by running the \code{estimate} function.} generic \code{forecast}.} \item{conditional_forecast}{a list of length \code{C} containing -\code{horizon x N} matrices or \code{ts} objects with of the same dimensions -with forecasted values for selected variables. These objects should only -contain \code{numeric} or \code{NA} values. The entries with \code{NA} values -correspond to the values that are forecasted conditionally on the realisations -provided as \code{numeric} values.} +\code{horizon x N} matrices with forecasted values for selected variables. +These matrices should only contain \code{numeric} or \code{NA} values. The +entries with \code{NA} values correspond to the values that are forecasted +conditionally on the realisations provided as \code{numeric} values.} } \value{ A list of class \code{PanelForecasts} containing the @@ -59,6 +58,25 @@ ilo_cubic_panel |> estimate(S = 20) |> forecast(horizon = 2) -> predictive +# conditional forecasting 6 years ahead conditioning on +# provided future values for the Gross Domestic Product +# growth rate +############################################################ +#' data(ilo_conditional_forecast) # load the conditional forecasts of dgdp +predictive = forecast(posterior, 6, conditional_forecast = ilo_conditional_forecast) + +# workflow with the pipe |> +############################################################ +set.seed(123) +ilo_cubic_panel |> + specify_bvarPANEL$new() |> + estimate(S = 10) |> + estimate(S = 20) |> + forecast( + horizon = 6, + conditional_forecast = ilo_conditional_forecast + ) -> predictive + } \author{ Tomasz Woźniak \email{wozniak.tom@pm.me} From 10c92e75d76f982e2572895f9394ddb60a426526 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tomasz=20Wo=C5=BAniak?= Date: Mon, 10 Jun 2024 23:07:11 +1000 Subject: [PATCH 8/8] testing for conditional forecasts #8 --- inst/tinytest/test_forecast.R | 43 +++++++++++++++++++++++++++++++++++ 1 file changed, 43 insertions(+) diff --git a/inst/tinytest/test_forecast.R b/inst/tinytest/test_forecast.R index 19f14bd..2d2750b 100644 --- a/inst/tinytest/test_forecast.R +++ b/inst/tinytest/test_forecast.R @@ -38,3 +38,46 @@ expect_error( forecast(run_no1, horizon = 1.5), info = "forecast: specify horizon as integer." ) + + +# conditional forecasting +data(ilo_conditional_forecast) + +set.seed(1) +suppressMessages( + specification_no1 <- specify_bvarPANEL$new(ilo_cubic_panel) +) +run_no1 <- estimate(specification_no1, 3, 1, show_progress = FALSE) +ff <- forecast(run_no1, 6, conditional_forecast = ilo_conditional_forecast) + +set.seed(1) +suppressMessages( + ff2 <- ilo_cubic_panel |> + specify_bvarPANEL$new() |> + estimate(S = 3, thin = 1, show_progress = FALSE) |> + forecast(horizon = 6, conditional_forecast = ilo_conditional_forecast) +) + + +expect_identical( + ff$forecasts[1,1,1,1], ff2$forecasts[1,1,1,1], + info = "conditional forecast: forecast identical for normal and pipe workflow." +) + +expect_true( + is.numeric(ff$forecasts) & is.array(ff$forecasts), + info = "conditional forecast: returns numeric array." +) + +expect_error( + forecast(run_no1, horizon = 4, conditional_forecast = ilo_conditional_forecast), + pattern = "horizon", + info = "conditional forecast: provided forecasts different from horizon." +) + +expect_error( + forecast(run_no1, horizon = 6, conditional_forecast = ilo_conditional_forecast[-1]), + info = "conditional forecast: uneven number of countries in forecasts and data." +) + +