diff --git a/DESCRIPTION b/DESCRIPTION new file mode 100644 index 0000000..a9c4e6a --- /dev/null +++ b/DESCRIPTION @@ -0,0 +1,38 @@ +Package: sovereign +Title: State-Dependant Empirical Analysis +Version: 0.0.0.9000 +Authors@R: + person(given = "Tyler", + family = "Pike", + role = c("aut", "cre"), + email = "tjpike7@gmail.com") +Description: A set of tools for state-dependant empirical analysis + through both VAR- and local projection-based state-dependant forecasts, + impulse response functions, and forecast error variance decomposition. +License: GPL-3 +Encoding: UTF-8 +LazyData: true +Roxygen: list(markdown = TRUE) +RoxygenNote: 7.1.1 +Imports: + broom, + dplyr, + ggplot2, + gridExtra, + lmtest, + lubridate, + magrittr, + purrr, + randomForest, + sandwich, + stats, + tidyr, + xts, + zoo +Suggests: + testthat, + knitr, + rmarkdown, + quantmod, + covr +VignetteBuilder: knitr diff --git a/LICENSE.md b/LICENSE.md new file mode 100644 index 0000000..e136d50 --- /dev/null +++ b/LICENSE.md @@ -0,0 +1,595 @@ +GNU General Public License +========================== + +_Version 3, 29 June 2007_ +_Copyright © 2007 Free Software Foundation, Inc. <>_ + +Everyone is permitted to copy and distribute verbatim copies of this license +document, but changing it is not allowed. + +## Preamble + +The GNU General Public License is a free, copyleft license for software and other +kinds of works. + +The licenses for most software and other practical works are designed to take away +your freedom to share and change the works. By contrast, the GNU General Public +License is intended to guarantee your freedom to share and change all versions of a +program--to make sure it remains free software for all its users. We, the Free +Software Foundation, use the GNU General Public License for most of our software; it +applies also to any other work released this way by its authors. You can apply it to +your programs, too. + +When we speak of free software, we are referring to freedom, not price. Our General +Public Licenses are designed to make sure that you have the freedom to distribute +copies of free software (and charge for them if you wish), that you receive source +code or can get it if you want it, that you can change the software or use pieces of +it in new free programs, and that you know you can do these things. + +To protect your rights, we need to prevent others from denying you these rights or +asking you to surrender the rights. Therefore, you have certain responsibilities if +you distribute copies of the software, or if you modify it: responsibilities to +respect the freedom of others. + +For example, if you distribute copies of such a program, whether gratis or for a fee, +you must pass on to the recipients the same freedoms that you received. You must make +sure that they, too, receive or can get the source code. And you must show them these +terms so they know their rights. + +Developers that use the GNU GPL protect your rights with two steps: **(1)** assert +copyright on the software, and **(2)** offer you this License giving you legal permission +to copy, distribute and/or modify it. + +For the developers' and authors' protection, the GPL clearly explains that there is +no warranty for this free software. For both users' and authors' sake, the GPL +requires that modified versions be marked as changed, so that their problems will not +be attributed erroneously to authors of previous versions. + +Some devices are designed to deny users access to install or run modified versions of +the software inside them, although the manufacturer can do so. This is fundamentally +incompatible with the aim of protecting users' freedom to change the software. The +systematic pattern of such abuse occurs in the area of products for individuals to +use, which is precisely where it is most unacceptable. Therefore, we have designed +this version of the GPL to prohibit the practice for those products. If such problems +arise substantially in other domains, we stand ready to extend this provision to +those domains in future versions of the GPL, as needed to protect the freedom of +users. + +Finally, every program is threatened constantly by software patents. States should +not allow patents to restrict development and use of software on general-purpose +computers, but in those that do, we wish to avoid the special danger that patents +applied to a free program could make it effectively proprietary. To prevent this, the +GPL assures that patents cannot be used to render the program non-free. + +The precise terms and conditions for copying, distribution and modification follow. + +## TERMS AND CONDITIONS + +### 0. Definitions + +“This License” refers to version 3 of the GNU General Public License. + +“Copyright” also means copyright-like laws that apply to other kinds of +works, such as semiconductor masks. + +“The Program” refers to any copyrightable work licensed under this +License. Each licensee is addressed as “you”. “Licensees” and +“recipients” may be individuals or organizations. + +To “modify” a work means to copy from or adapt all or part of the work in +a fashion requiring copyright permission, other than the making of an exact copy. The +resulting work is called a “modified version” of the earlier work or a +work “based on” the earlier work. + +A “covered work” means either the unmodified Program or a work based on +the Program. + +To “propagate” a work means to do anything with it that, without +permission, would make you directly or secondarily liable for infringement under +applicable copyright law, except executing it on a computer or modifying a private +copy. Propagation includes copying, distribution (with or without modification), +making available to the public, and in some countries other activities as well. + +To “convey” a work means any kind of propagation that enables other +parties to make or receive copies. Mere interaction with a user through a computer +network, with no transfer of a copy, is not conveying. + +An interactive user interface displays “Appropriate Legal Notices” to the +extent that it includes a convenient and prominently visible feature that **(1)** +displays an appropriate copyright notice, and **(2)** tells the user that there is no +warranty for the work (except to the extent that warranties are provided), that +licensees may convey the work under this License, and how to view a copy of this +License. If the interface presents a list of user commands or options, such as a +menu, a prominent item in the list meets this criterion. + +### 1. Source Code + +The “source code” for a work means the preferred form of the work for +making modifications to it. “Object code” means any non-source form of a +work. + +A “Standard Interface” means an interface that either is an official +standard defined by a recognized standards body, or, in the case of interfaces +specified for a particular programming language, one that is widely used among +developers working in that language. + +The “System Libraries” of an executable work include anything, other than +the work as a whole, that **(a)** is included in the normal form of packaging a Major +Component, but which is not part of that Major Component, and **(b)** serves only to +enable use of the work with that Major Component, or to implement a Standard +Interface for which an implementation is available to the public in source code form. +A “Major Component”, in this context, means a major essential component +(kernel, window system, and so on) of the specific operating system (if any) on which +the executable work runs, or a compiler used to produce the work, or an object code +interpreter used to run it. + +The “Corresponding Source” for a work in object code form means all the +source code needed to generate, install, and (for an executable work) run the object +code and to modify the work, including scripts to control those activities. However, +it does not include the work's System Libraries, or general-purpose tools or +generally available free programs which are used unmodified in performing those +activities but which are not part of the work. For example, Corresponding Source +includes interface definition files associated with source files for the work, and +the source code for shared libraries and dynamically linked subprograms that the work +is specifically designed to require, such as by intimate data communication or +control flow between those subprograms and other parts of the work. + +The Corresponding Source need not include anything that users can regenerate +automatically from other parts of the Corresponding Source. + +The Corresponding Source for a work in source code form is that same work. + +### 2. Basic Permissions + +All rights granted under this License are granted for the term of copyright on the +Program, and are irrevocable provided the stated conditions are met. This License +explicitly affirms your unlimited permission to run the unmodified Program. The +output from running a covered work is covered by this License only if the output, +given its content, constitutes a covered work. This License acknowledges your rights +of fair use or other equivalent, as provided by copyright law. + +You may make, run and propagate covered works that you do not convey, without +conditions so long as your license otherwise remains in force. You may convey covered +works to others for the sole purpose of having them make modifications exclusively +for you, or provide you with facilities for running those works, provided that you +comply with the terms of this License in conveying all material for which you do not +control copyright. Those thus making or running the covered works for you must do so +exclusively on your behalf, under your direction and control, on terms that prohibit +them from making any copies of your copyrighted material outside their relationship +with you. + +Conveying under any other circumstances is permitted solely under the conditions +stated below. Sublicensing is not allowed; section 10 makes it unnecessary. + +### 3. Protecting Users' Legal Rights From Anti-Circumvention Law + +No covered work shall be deemed part of an effective technological measure under any +applicable law fulfilling obligations under article 11 of the WIPO copyright treaty +adopted on 20 December 1996, or similar laws prohibiting or restricting circumvention +of such measures. + +When you convey a covered work, you waive any legal power to forbid circumvention of +technological measures to the extent such circumvention is effected by exercising +rights under this License with respect to the covered work, and you disclaim any +intention to limit operation or modification of the work as a means of enforcing, +against the work's users, your or third parties' legal rights to forbid circumvention +of technological measures. + +### 4. Conveying Verbatim Copies + +You may convey verbatim copies of the Program's source code as you receive it, in any +medium, provided that you conspicuously and appropriately publish on each copy an +appropriate copyright notice; keep intact all notices stating that this License and +any non-permissive terms added in accord with section 7 apply to the code; keep +intact all notices of the absence of any warranty; and give all recipients a copy of +this License along with the Program. + +You may charge any price or no price for each copy that you convey, and you may offer +support or warranty protection for a fee. + +### 5. Conveying Modified Source Versions + +You may convey a work based on the Program, or the modifications to produce it from +the Program, in the form of source code under the terms of section 4, provided that +you also meet all of these conditions: + +* **a)** The work must carry prominent notices stating that you modified it, and giving a +relevant date. +* **b)** The work must carry prominent notices stating that it is released under this +License and any conditions added under section 7. This requirement modifies the +requirement in section 4 to “keep intact all notices”. +* **c)** You must license the entire work, as a whole, under this License to anyone who +comes into possession of a copy. This License will therefore apply, along with any +applicable section 7 additional terms, to the whole of the work, and all its parts, +regardless of how they are packaged. This License gives no permission to license the +work in any other way, but it does not invalidate such permission if you have +separately received it. +* **d)** If the work has interactive user interfaces, each must display Appropriate Legal +Notices; however, if the Program has interactive interfaces that do not display +Appropriate Legal Notices, your work need not make them do so. + +A compilation of a covered work with other separate and independent works, which are +not by their nature extensions of the covered work, and which are not combined with +it such as to form a larger program, in or on a volume of a storage or distribution +medium, is called an “aggregate” if the compilation and its resulting +copyright are not used to limit the access or legal rights of the compilation's users +beyond what the individual works permit. Inclusion of a covered work in an aggregate +does not cause this License to apply to the other parts of the aggregate. + +### 6. Conveying Non-Source Forms + +You may convey a covered work in object code form under the terms of sections 4 and +5, provided that you also convey the machine-readable Corresponding Source under the +terms of this License, in one of these ways: + +* **a)** Convey the object code in, or embodied in, a physical product (including a +physical distribution medium), accompanied by the Corresponding Source fixed on a +durable physical medium customarily used for software interchange. +* **b)** Convey the object code in, or embodied in, a physical product (including a +physical distribution medium), accompanied by a written offer, valid for at least +three years and valid for as long as you offer spare parts or customer support for +that product model, to give anyone who possesses the object code either **(1)** a copy of +the Corresponding Source for all the software in the product that is covered by this +License, on a durable physical medium customarily used for software interchange, for +a price no more than your reasonable cost of physically performing this conveying of +source, or **(2)** access to copy the Corresponding Source from a network server at no +charge. +* **c)** Convey individual copies of the object code with a copy of the written offer to +provide the Corresponding Source. This alternative is allowed only occasionally and +noncommercially, and only if you received the object code with such an offer, in +accord with subsection 6b. +* **d)** Convey the object code by offering access from a designated place (gratis or for +a charge), and offer equivalent access to the Corresponding Source in the same way +through the same place at no further charge. You need not require recipients to copy +the Corresponding Source along with the object code. If the place to copy the object +code is a network server, the Corresponding Source may be on a different server +(operated by you or a third party) that supports equivalent copying facilities, +provided you maintain clear directions next to the object code saying where to find +the Corresponding Source. Regardless of what server hosts the Corresponding Source, +you remain obligated to ensure that it is available for as long as needed to satisfy +these requirements. +* **e)** Convey the object code using peer-to-peer transmission, provided you inform +other peers where the object code and Corresponding Source of the work are being +offered to the general public at no charge under subsection 6d. + +A separable portion of the object code, whose source code is excluded from the +Corresponding Source as a System Library, need not be included in conveying the +object code work. + +A “User Product” is either **(1)** a “consumer product”, which +means any tangible personal property which is normally used for personal, family, or +household purposes, or **(2)** anything designed or sold for incorporation into a +dwelling. In determining whether a product is a consumer product, doubtful cases +shall be resolved in favor of coverage. For a particular product received by a +particular user, “normally used” refers to a typical or common use of +that class of product, regardless of the status of the particular user or of the way +in which the particular user actually uses, or expects or is expected to use, the +product. A product is a consumer product regardless of whether the product has +substantial commercial, industrial or non-consumer uses, unless such uses represent +the only significant mode of use of the product. + +“Installation Information” for a User Product means any methods, +procedures, authorization keys, or other information required to install and execute +modified versions of a covered work in that User Product from a modified version of +its Corresponding Source. The information must suffice to ensure that the continued +functioning of the modified object code is in no case prevented or interfered with +solely because modification has been made. + +If you convey an object code work under this section in, or with, or specifically for +use in, a User Product, and the conveying occurs as part of a transaction in which +the right of possession and use of the User Product is transferred to the recipient +in perpetuity or for a fixed term (regardless of how the transaction is +characterized), the Corresponding Source conveyed under this section must be +accompanied by the Installation Information. But this requirement does not apply if +neither you nor any third party retains the ability to install modified object code +on the User Product (for example, the work has been installed in ROM). + +The requirement to provide Installation Information does not include a requirement to +continue to provide support service, warranty, or updates for a work that has been +modified or installed by the recipient, or for the User Product in which it has been +modified or installed. Access to a network may be denied when the modification itself +materially and adversely affects the operation of the network or violates the rules +and protocols for communication across the network. + +Corresponding Source conveyed, and Installation Information provided, in accord with +this section must be in a format that is publicly documented (and with an +implementation available to the public in source code form), and must require no +special password or key for unpacking, reading or copying. + +### 7. Additional Terms + +“Additional permissions” are terms that supplement the terms of this +License by making exceptions from one or more of its conditions. Additional +permissions that are applicable to the entire Program shall be treated as though they +were included in this License, to the extent that they are valid under applicable +law. If additional permissions apply only to part of the Program, that part may be +used separately under those permissions, but the entire Program remains governed by +this License without regard to the additional permissions. + +When you convey a copy of a covered work, you may at your option remove any +additional permissions from that copy, or from any part of it. (Additional +permissions may be written to require their own removal in certain cases when you +modify the work.) You may place additional permissions on material, added by you to a +covered work, for which you have or can give appropriate copyright permission. + +Notwithstanding any other provision of this License, for material you add to a +covered work, you may (if authorized by the copyright holders of that material) +supplement the terms of this License with terms: + +* **a)** Disclaiming warranty or limiting liability differently from the terms of +sections 15 and 16 of this License; or +* **b)** Requiring preservation of specified reasonable legal notices or author +attributions in that material or in the Appropriate Legal Notices displayed by works +containing it; or +* **c)** Prohibiting misrepresentation of the origin of that material, or requiring that +modified versions of such material be marked in reasonable ways as different from the +original version; or +* **d)** Limiting the use for publicity purposes of names of licensors or authors of the +material; or +* **e)** Declining to grant rights under trademark law for use of some trade names, +trademarks, or service marks; or +* **f)** Requiring indemnification of licensors and authors of that material by anyone +who conveys the material (or modified versions of it) with contractual assumptions of +liability to the recipient, for any liability that these contractual assumptions +directly impose on those licensors and authors. + +All other non-permissive additional terms are considered “further +restrictions” within the meaning of section 10. If the Program as you received +it, or any part of it, contains a notice stating that it is governed by this License +along with a term that is a further restriction, you may remove that term. If a +license document contains a further restriction but permits relicensing or conveying +under this License, you may add to a covered work material governed by the terms of +that license document, provided that the further restriction does not survive such +relicensing or conveying. + +If you add terms to a covered work in accord with this section, you must place, in +the relevant source files, a statement of the additional terms that apply to those +files, or a notice indicating where to find the applicable terms. + +Additional terms, permissive or non-permissive, may be stated in the form of a +separately written license, or stated as exceptions; the above requirements apply +either way. + +### 8. Termination + +You may not propagate or modify a covered work except as expressly provided under +this License. Any attempt otherwise to propagate or modify it is void, and will +automatically terminate your rights under this License (including any patent licenses +granted under the third paragraph of section 11). + +However, if you cease all violation of this License, then your license from a +particular copyright holder is reinstated **(a)** provisionally, unless and until the +copyright holder explicitly and finally terminates your license, and **(b)** permanently, +if the copyright holder fails to notify you of the violation by some reasonable means +prior to 60 days after the cessation. + +Moreover, your license from a particular copyright holder is reinstated permanently +if the copyright holder notifies you of the violation by some reasonable means, this +is the first time you have received notice of violation of this License (for any +work) from that copyright holder, and you cure the violation prior to 30 days after +your receipt of the notice. + +Termination of your rights under this section does not terminate the licenses of +parties who have received copies or rights from you under this License. If your +rights have been terminated and not permanently reinstated, you do not qualify to +receive new licenses for the same material under section 10. + +### 9. Acceptance Not Required for Having Copies + +You are not required to accept this License in order to receive or run a copy of the +Program. Ancillary propagation of a covered work occurring solely as a consequence of +using peer-to-peer transmission to receive a copy likewise does not require +acceptance. However, nothing other than this License grants you permission to +propagate or modify any covered work. These actions infringe copyright if you do not +accept this License. Therefore, by modifying or propagating a covered work, you +indicate your acceptance of this License to do so. + +### 10. Automatic Licensing of Downstream Recipients + +Each time you convey a covered work, the recipient automatically receives a license +from the original licensors, to run, modify and propagate that work, subject to this +License. You are not responsible for enforcing compliance by third parties with this +License. + +An “entity transaction” is a transaction transferring control of an +organization, or substantially all assets of one, or subdividing an organization, or +merging organizations. If propagation of a covered work results from an entity +transaction, each party to that transaction who receives a copy of the work also +receives whatever licenses to the work the party's predecessor in interest had or +could give under the previous paragraph, plus a right to possession of the +Corresponding Source of the work from the predecessor in interest, if the predecessor +has it or can get it with reasonable efforts. + +You may not impose any further restrictions on the exercise of the rights granted or +affirmed under this License. For example, you may not impose a license fee, royalty, +or other charge for exercise of rights granted under this License, and you may not +initiate litigation (including a cross-claim or counterclaim in a lawsuit) alleging +that any patent claim is infringed by making, using, selling, offering for sale, or +importing the Program or any portion of it. + +### 11. Patents + +A “contributor” is a copyright holder who authorizes use under this +License of the Program or a work on which the Program is based. The work thus +licensed is called the contributor's “contributor version”. + +A contributor's “essential patent claims” are all patent claims owned or +controlled by the contributor, whether already acquired or hereafter acquired, that +would be infringed by some manner, permitted by this License, of making, using, or +selling its contributor version, but do not include claims that would be infringed +only as a consequence of further modification of the contributor version. For +purposes of this definition, “control” includes the right to grant patent +sublicenses in a manner consistent with the requirements of this License. + +Each contributor grants you a non-exclusive, worldwide, royalty-free patent license +under the contributor's essential patent claims, to make, use, sell, offer for sale, +import and otherwise run, modify and propagate the contents of its contributor +version. + +In the following three paragraphs, a “patent license” is any express +agreement or commitment, however denominated, not to enforce a patent (such as an +express permission to practice a patent or covenant not to sue for patent +infringement). To “grant” such a patent license to a party means to make +such an agreement or commitment not to enforce a patent against the party. + +If you convey a covered work, knowingly relying on a patent license, and the +Corresponding Source of the work is not available for anyone to copy, free of charge +and under the terms of this License, through a publicly available network server or +other readily accessible means, then you must either **(1)** cause the Corresponding +Source to be so available, or **(2)** arrange to deprive yourself of the benefit of the +patent license for this particular work, or **(3)** arrange, in a manner consistent with +the requirements of this License, to extend the patent license to downstream +recipients. “Knowingly relying” means you have actual knowledge that, but +for the patent license, your conveying the covered work in a country, or your +recipient's use of the covered work in a country, would infringe one or more +identifiable patents in that country that you have reason to believe are valid. + +If, pursuant to or in connection with a single transaction or arrangement, you +convey, or propagate by procuring conveyance of, a covered work, and grant a patent +license to some of the parties receiving the covered work authorizing them to use, +propagate, modify or convey a specific copy of the covered work, then the patent +license you grant is automatically extended to all recipients of the covered work and +works based on it. + +A patent license is “discriminatory” if it does not include within the +scope of its coverage, prohibits the exercise of, or is conditioned on the +non-exercise of one or more of the rights that are specifically granted under this +License. You may not convey a covered work if you are a party to an arrangement with +a third party that is in the business of distributing software, under which you make +payment to the third party based on the extent of your activity of conveying the +work, and under which the third party grants, to any of the parties who would receive +the covered work from you, a discriminatory patent license **(a)** in connection with +copies of the covered work conveyed by you (or copies made from those copies), or **(b)** +primarily for and in connection with specific products or compilations that contain +the covered work, unless you entered into that arrangement, or that patent license +was granted, prior to 28 March 2007. + +Nothing in this License shall be construed as excluding or limiting any implied +license or other defenses to infringement that may otherwise be available to you +under applicable patent law. + +### 12. No Surrender of Others' Freedom + +If conditions are imposed on you (whether by court order, agreement or otherwise) +that contradict the conditions of this License, they do not excuse you from the +conditions of this License. If you cannot convey a covered work so as to satisfy +simultaneously your obligations under this License and any other pertinent +obligations, then as a consequence you may not convey it at all. For example, if you +agree to terms that obligate you to collect a royalty for further conveying from +those to whom you convey the Program, the only way you could satisfy both those terms +and this License would be to refrain entirely from conveying the Program. + +### 13. Use with the GNU Affero General Public License + +Notwithstanding any other provision of this License, you have permission to link or +combine any covered work with a work licensed under version 3 of the GNU Affero +General Public License into a single combined work, and to convey the resulting work. +The terms of this License will continue to apply to the part which is the covered +work, but the special requirements of the GNU Affero General Public License, section +13, concerning interaction through a network will apply to the combination as such. + +### 14. Revised Versions of this License + +The Free Software Foundation may publish revised and/or new versions of the GNU +General Public License from time to time. Such new versions will be similar in spirit +to the present version, but may differ in detail to address new problems or concerns. + +Each version is given a distinguishing version number. If the Program specifies that +a certain numbered version of the GNU General Public License “or any later +version” applies to it, you have the option of following the terms and +conditions either of that numbered version or of any later version published by the +Free Software Foundation. If the Program does not specify a version number of the GNU +General Public License, you may choose any version ever published by the Free +Software Foundation. + +If the Program specifies that a proxy can decide which future versions of the GNU +General Public License can be used, that proxy's public statement of acceptance of a +version permanently authorizes you to choose that version for the Program. + +Later license versions may give you additional or different permissions. However, no +additional obligations are imposed on any author or copyright holder as a result of +your choosing to follow a later version. + +### 15. Disclaimer of Warranty + +THERE IS NO WARRANTY FOR THE PROGRAM, TO THE EXTENT PERMITTED BY APPLICABLE LAW. +EXCEPT WHEN OTHERWISE STATED IN WRITING THE COPYRIGHT HOLDERS AND/OR OTHER PARTIES +PROVIDE THE PROGRAM “AS IS” WITHOUT WARRANTY OF ANY KIND, EITHER +EXPRESSED OR IMPLIED, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF +MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. THE ENTIRE RISK AS TO THE +QUALITY AND PERFORMANCE OF THE PROGRAM IS WITH YOU. SHOULD THE PROGRAM PROVE +DEFECTIVE, YOU ASSUME THE COST OF ALL NECESSARY SERVICING, REPAIR OR CORRECTION. + +### 16. Limitation of Liability + +IN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW OR AGREED TO IN WRITING WILL ANY +COPYRIGHT HOLDER, OR ANY OTHER PARTY WHO MODIFIES AND/OR CONVEYS THE PROGRAM AS +PERMITTED ABOVE, BE LIABLE TO YOU FOR DAMAGES, INCLUDING ANY GENERAL, SPECIAL, +INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING OUT OF THE USE OR INABILITY TO USE THE +PROGRAM (INCLUDING BUT NOT LIMITED TO LOSS OF DATA OR DATA BEING RENDERED INACCURATE +OR LOSSES SUSTAINED BY YOU OR THIRD PARTIES OR A FAILURE OF THE PROGRAM TO OPERATE +WITH ANY OTHER PROGRAMS), EVEN IF SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF THE +POSSIBILITY OF SUCH DAMAGES. + +### 17. Interpretation of Sections 15 and 16 + +If the disclaimer of warranty and limitation of liability provided above cannot be +given local legal effect according to their terms, reviewing courts shall apply local +law that most closely approximates an absolute waiver of all civil liability in +connection with the Program, unless a warranty or assumption of liability accompanies +a copy of the Program in return for a fee. + +_END OF TERMS AND CONDITIONS_ + +## How to Apply These Terms to Your New Programs + +If you develop a new program, and you want it to be of the greatest possible use to +the public, the best way to achieve this is to make it free software which everyone +can redistribute and change under these terms. + +To do so, attach the following notices to the program. It is safest to attach them +to the start of each source file to most effectively state the exclusion of warranty; +and each file should have at least the “copyright” line and a pointer to +where the full notice is found. + + + Copyright (C) 2021 Tyler J. Pike + + This program is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with this program. If not, see . + +Also add information on how to contact you by electronic and paper mail. + +If the program does terminal interaction, make it output a short notice like this +when it starts in an interactive mode: + + sovereign Copyright (C) 2021 Tyler J. Pike + This program comes with ABSOLUTELY NO WARRANTY; for details type 'show w'. + This is free software, and you are welcome to redistribute it + under certain conditions; type 'show c' for details. + +The hypothetical commands `show w` and `show c` should show the appropriate parts of +the General Public License. Of course, your program's commands might be different; +for a GUI interface, you would use an “about box”. + +You should also get your employer (if you work as a programmer) or school, if any, to +sign a “copyright disclaimer” for the program, if necessary. For more +information on this, and how to apply and follow the GNU GPL, see +<>. + +The GNU General Public License does not permit incorporating your program into +proprietary programs. If your program is a subroutine library, you may consider it +more useful to permit linking proprietary applications with the library. If this is +what you want to do, use the GNU Lesser General Public License instead of this +License. But first, please read +<>. diff --git a/NAMESPACE b/NAMESPACE new file mode 100644 index 0000000..3d375d5 --- /dev/null +++ b/NAMESPACE @@ -0,0 +1,16 @@ +# Generated by roxygen2: do not edit by hand + +export("%>%") +export(VAR) +export(individual_var_irf_plot) +export(learn_regimes) +export(lp_irf) +export(lp_irf_chart) +export(threshold_VAR) +export(threshold_lp_irf) +export(threshold_var_fevd) +export(threshold_var_irf) +export(var_fevd) +export(var_irf) +export(var_irf_plot) +importFrom(magrittr,"%>%") diff --git a/R/external_imports.R b/R/external_imports.R new file mode 100644 index 0000000..e79f3d8 --- /dev/null +++ b/R/external_imports.R @@ -0,0 +1,11 @@ +#' Pipe operator +#' +#' See \code{magrittr::\link[magrittr:pipe]{\%>\%}} for details. +#' +#' @name %>% +#' @rdname pipe +#' @keywords internal +#' @export +#' @importFrom magrittr %>% +#' @usage lhs \%>\% rhs +NULL diff --git a/R/helper.R b/R/helper.R new file mode 100644 index 0000000..4562e65 --- /dev/null +++ b/R/helper.R @@ -0,0 +1,58 @@ +#------------------------------------------ +# Data helper functions +# (copied from OOS package) +#------------------------------------------ +# create n lags +n.lag = function( + Data, # data.frame: data frame of variables to lag and a 'date' column + lags, # int: number of lags to create + variables = NULL # string: vector of variable names to lag, default is all non-date variables +){ + + if(is.null(variables)){ + variables = names(dplyr::select(Data, -contains('date'))) + } + + Data = c(0:lags) %>% + purrr::map( + .f = function(n){ + + if(n == 0){return(Data)} + + X = Data %>% + dplyr::mutate_at(variables, dplyr::lag, n) + + names(X)[names(X) != 'date'] = paste0(names(X)[names(X) != 'date'], '.l', n) + + return(X) + } + ) %>% + purrr::reduce(dplyr::full_join, by = 'date') + + + return(Data) +} + +# adjust forecast dates +forecast_date = function( + forecast.date, + horizon, + freq +){ + + date = forecast.date + + if(freq == 'day'){ + date = forecast.date + horizon + }else if(freq == 'week'){ + lubridate::week(date) = lubridate::week(date) + horizon + }else if(freq == 'month'){ + lubridate::month(date) = lubridate::month(date) + horizon + }else if(freq == 'quarter'){ + lubridate::month(date) = lubridate::month(date) + horizon*3 + }else if(freq == 'year'){ + lubridate::year(date) = lubridate::year(date) + horizon + } + + return(date) +} \ No newline at end of file diff --git a/R/local_projections.R b/R/local_projections.R new file mode 100644 index 0000000..156cc5a --- /dev/null +++ b/R/local_projections.R @@ -0,0 +1,242 @@ + +#------------------------------------------------- +# Function to produce IRF +#------------------------------------------------- +#' Estimate single-regime local projection IRFs +#' +#' @param data data.frame, matrix, ts, xts, zoo: Endogenous regressors +#' @param shock string: variable to shock +#' @param target string: variable betas to collect +#' @param horizons int: horizons to forecast out to +#' @param lags int: lags to include in regressions +#' +#' @return data.frame +#' +#' @examples +#' \dontrun{ +#' lp_irf( +#' data = Data, +#' shock = 'x', +#' target = 'y', +#' horizons = 20, +#' lags = 2) +#' } +#' +#' @export + +lp_irf = function( + data, # dataframe of covariate + shock, # string denoting variable to shock + target, # string denoting variable betas to collect + horizons, # horizons to forecast out to + lags # lags to include in regressions +){ + + # function warnings + if(!is.matrix(data) & !is.data.frame(data)){ + errorCondition('data must be a matrix or data.frame') + } + if(!is.numeric(lags) | lags %% 1 != 0){ + errorCondition('lags must be an integer') + } + if(!is.numeric(horizons) | horizons %% 1 != 0 | horizons <= 0){ + errorCondition('horizons must be a positive integer') + } + + # cast as data frame if ts, xts, or zoo object + if(is.ts(data) | xts::is.xts(data) | zoo::is.zoo(data)){ + data = data.frame(date = zoo::index(date), data) + } + + # first create the proper table and variable names + data = data %>% dplyr::rename(target = target) + data = as.data.frame(data) + data = na.omit(data) + + #################################################### + # generate the lags and leads + #################################################### + date = data$date + data = data %>% dplyr::select(-date) + final = data + # generate lags + for(i in 1:lags){ + temp = as.data.frame(lapply(data, MARGIN = 2, FUN = lag, n=i)) + colnames(temp) = paste0(colnames(temp),'.lag',i) + final = cbind(temp,final) + } + # generate leads + # not the most efficient but it works nonetheless + temp = matrix(ncol = horizons, nrow = length(data$target)) + for(i in 1:horizons){ + temp[,i] = dplyr::lead(data$target, n = i) + } + colnames(temp) = paste0('target.l',c(1:horizons)) + leads = colnames(temp) #used later + final = cbind(temp,final) + final$date = date + data = final + + + #################################################### + # perform local impulse response regressions + #################################################### + # strip out the non-explanetory variables + explanetory = data %>% dplyr::select(-date, -leads) + + # strip out target variables + targets = data %>% dplyr::select(leads) + + # generate the coef and sd for plotting + # create storage matrix + irfData = matrix(ncol = 3, nrow = horizons) + colnames(irfData) = c('Horizon','Coef','Std.dev') + irfData[,1] = c(1:horizons) + # calculate regressions + for(i in 1:horizons){ + # generate the projections + directProjection = stats::lm(targets[,i] ~., data = explanetory) + out = lmtest::coeftest(directProjection,vcov = sandwich::NeweyWest(directProjection,prewhite=FALSE)) + shockIndex = match(shock,rownames(out)) + # store the data + irfData[i,2] <- out[shockIndex,1] + irfData[i,3] <- out[shockIndex,2] + } + + #################################################### + # finalize data and return + #################################################### + irfData = as.data.frame(irfData) + # generate and store upper/lower bound + irfData$lowerBound <- irfData$Coef - 1.64*irfData$Std + irfData$upperBound <- irfData$Coef + 1.64*irfData$Std + return(irfData) + +} + +#------------------------------------------------- +# Function to produce threshold IRF +# calculates responses in n different regimes +#------------------------------------------------- +#' Estimate multi-regime local projection IRFs +#' +#' @param data data.frame, matrix, ts, xts, zoo: Endogenous regressors +#' @param shock string: variable to shock +#' @param target string: variable betas to collect +#' @param thresholdVar data.frame of regime binaries +#' @param horizons int: horizons to forecast out to +#' @param lags int: lags to include in regressions +#' +#' @return data.frame +#' +#' @examples +#' \dontrun{ +#' threshold_lp_irf( +#' data = Data, +#' shock = 'x', +#' target = 'y', +#' thresholdVar = Data.regime, +#' horizons = 20, +#' lags = 2) +#' } +#' +#' @export + +threshold_lp_irf <- function( + data, # dataframe of covariate + shock, # string denoting variable to shock + target, # string denoting variable betas to collect + thresholdVar = NULL, # dataframe of regime binaries + horizons, # horizons to forecast out to + lags){ # lags to include in regressions + + # first create the proper table and variable names + data = data %>% dplyr::rename(target = target) + data = as.data.frame(data) + data = na.omit(data) + + # use only dates in both data and thresholdvar + data = data %>% dplyr::filter(date %in% thresholdVar$date) + thresholdVar = thresholdVar %>% dplyr::filter(date %in% data$date) + + #################################################### + # generate the lags and leads + #################################################### + date = data$date + data = data %>% dplyr::select(-date) + final = data + # generate lags + for(i in 1:lags){ + temp = as.data.frame(lapply(data, MARGIN = 2, FUN = lag, n=i)) + colnames(temp) = paste0(colnames(temp),'.lag',i) + final = cbind(temp,final) + } + # generate leads + # not the most efficient but it works nonetheless + temp = matrix(ncol = horizons, nrow = length(data$target)) + for(i in 1:horizons){ + temp[,i] = dplyr::lead(data$target, n = i) + } + colnames(temp) = paste0('target.l',c(1:horizons)) + leads = colnames(temp) #used later + final = cbind(temp,final) + final$date = date + data = final + + + #################################################### + # perform local impulse response regressions + #################################################### + # strip out the non-explanetory variables + explanetory = data %>% dplyr::select(-leads) + + # strip out target variables + targets = data %>% dplyr::select(leads) + + # create list to store results per threshold + outputList = list() + + # iteratre through thresholds creating the IRs + for(t in 1:(ncol(thresholdVar)-1)){ + # threshold variable + threshold = dplyr::select(thresholdVar, date, colnames(dplyr::select(thresholdVar, -date))[t]) %>% as.data.frame() + # carefully align s.t. matrix goes date, threshold, explanetory variables + X = dplyr::inner_join(threshold, explanetory, by = 'date') + # condition each row on regime (this appears to be the most robust way to do it) + X = data.frame(t(t(X[,3:ncol(X)] * X[,2]))) + + # generate the coef and sd for plotting + # create storage matrix + irfData = matrix(ncol = 3, nrow = horizons) + colnames(irfData) = c('Horizon','Coef','Std.dev') + irfData[,1] = c(1:horizons) + # calculate regressions + for(i in 1:horizons){ + # generate the projections + directProjection = lm(targets[,i] ~., data = X) + out = lmtest::coeftest(directProjection, + vcov = sandwich::NeweyWest(directProjection, + lag = lags, + prewhite=FALSE)) + shockIndex = match(shock,rownames(out)) + # store the data + irfData[i,2] = out[shockIndex,1] + irfData[i,3] = out[shockIndex,2] + } + + # finalize data and place in list + irfData = as.data.frame(irfData) + # generate and store upper/lower bound + irfData$lowerBound <- irfData$Coef - 1.64*irfData$Std + irfData$upperBound <- irfData$Coef + 1.64*irfData$Std + + outputList[[t]] = irfData + } + + #################################################### + # finalizing and returning output + #################################################### + names(outputList) = colnames(dplyr::select(thresholdVar, -date)) + return(outputList) +} + diff --git a/R/lp_irf_chart.R b/R/lp_irf_chart.R new file mode 100644 index 0000000..5d3be47 --- /dev/null +++ b/R/lp_irf_chart.R @@ -0,0 +1,223 @@ +# File: direct_projection_chart_functions.R +# Author: Tyler Pike +# Section: MA-MFA +# Date: 2/5/2020 +# Note(s): Houses functions to chart different types of Jorda (2005) style direct projections + + +#------------------------------------------------- +# Function to plot IRF +#------------------------------------------------- +#' Chart local projection IRFs +#' +#' @param plotData lp_irf output +#' @param title string: chart title +#' @param ylim int: y-axis limits, c(lower limit, upper limit) +#' @param ystep int: step size inbetween y-axis tick marks +#' @param ylab string: y-axis label +#' +#' @return plot +#' +#' @export + +lp_irf_chart = function( + plotData, + title, + ylim, + ystep, + ylab){ + + #create the functions and set options for creating the plots + opt = list(frame.lwd = 1.5, + line.lwd = c(1.7, 2), + label.cex = .75, + axis.cex = .65, + exhibit.title.cex = 1, + chart.title.line = 0.8, + chart.title.cex = 0.8, + foot.cex = .55, + legend.cex = .7, + line.types = rep(1,20), + legend.inset = .03, + yaxis.line = 0, + tick.length = 0.025, + yaxis.pos= .5, + colors = c("black","firebrick","SteelBlue", "DarkOliveGreen3", "goldenrod1", "blueviolet","magenta"), + key_adj_x=1, + key.cex=0.75, + paneltitleline=1.2, + paneltitlecex=.85, + tealbook.months = c("Jan.", "Feb.", "Mar.", "Apr.", "May", "Jun.", "July", "Aug.", "Sept.", "Oct.", "Nov.", "Dec."), + tealbook.quarters = c("Q1","Q2","Q3","Q4")) + + + + forecasts <- function(dbpath=NULL, + keynames="", + chart_title="", + units="", + ymin, + ymax, + ystep, + xtickfreq, + frequency="", + horizshift=0, + vertshift=0, + note="", + footnote.placement=1, + key="topleft", + start.date=min(plot.data$date), + end.date=max(plot.data$date), + show.recent.months=FALSE, + n.months=0, + mult.conv.opt=1, + legcol=1, + lineatzero=FALSE, + dbtype="fame", + dataframe=NULL, + extra.xlim = 0, + colors=opt$colors, + interp=TRUE, + xmajortickfreq="year", + xminortickfreq="month", + xhighfreq="month", + xlowfreq="year", + hadj=opt$yaxis.pos, + lty1=opt$line.types[1], + lty2=opt$line.types[1], + two.printlastvalue=FALSE, + horizshift2=0, + vertshift2=0, + rmlastlabel=FALSE, + lowdateplacement=seq(place1, place2, by = xlowfreq), + xmin=ymin, + xmax=ymax, + xstep=ystep, + yaxis.shock.label="") { + + xmin = 0 + xmax = nrow(plotData) + xstep = 1 + + + plot(plotData$Horizon, plotData$Coef, + type = 'l', + xlab = "", + ylab = "", + ylim = ylim,#c(ymin,ymax), + axes = FALSE, + col = "firebrick", + #lty = lty1, + #lwd = 1.7, + main = "", + bty = 'u', + yaxs = "i", + xaxs = "i", + xlim = c(.8,(nrow(plotData) + .2)), + pch = 16 #Consider 16, 19, 20 + #cex=1.3 + ) + + + polygon(c(plotData$Horizon,rev(plotData$Horizon)),c(plotData$upperBound,rev(plotData$lowerBound)),col="lightblue",border=NA) + + par(new=TRUE) + plot(plotData$Horizon, plotData$Coef, + type = 'l', + xlab = "", + ylab = "", + ylim = ylim,#c(ymin,ymax), + axes = FALSE, + col = "firebrick", + #lty = lty1, + #lwd = 1.7, + main = "", + bty = 'u', + yaxs = "i", + xaxs = "i", + xlim = c(.8,(nrow(plotData) + .2)), + pch = 16 #Consider 16, 19, 20 + #cex=1.3 + ) + + + set_chart_parameters() + plotHookBox(lwd=opt$frame.lwd) + + #y axis + + axis(side = 4, + at = seq(ymin,ymax, by = ystep), + tck = opt$tick.length, + cex.axis = opt$axis.cex, + las = 2, + hadj=hadj) + + axis(side = 2, + at = seq(ymin,ymax, by = ystep), + tck = opt$tick.length, + cex.axis = opt$axis.cex, + las = 2, + labels=FALSE) + + + axis(side = 1, at = seq(xmin,xmax,by = xstep), + tck = opt$tick.length+.015, + cex.axis = opt$axis.cex, + las = 0, + labels=FALSE,hadj=hadj) + + + axis(side = 1, at = seq(xmin,xmax,by = xstep), + tick=FALSE, + cex.axis = opt$axis.cex, + las = 0, + labels=TRUE, + hadj=hadj, + line = -1) + + + #title + # + # title(main = chart_title, + # cex.main = paneltitle.cex, + # line=.1,font.main=1, + # adj=0) + + #labels + mtext(units, side = 3, line =0 , adj = 1, outer = FALSE, cex = opt$legend.cex) + mtext(chart_title, side = 3, line =0 , adj = 0, outer = FALSE, cex = paneltitle.cex) + mtext(frequency, side = 3, line =-.5 , adj = .08, outer = FALSE, cex = opt$legend.cex) + mtext(note, side = 1, line = footnote.placement, adj = 0, outer = FALSE, cex = opt$legend.cex) + mtext(yaxis.shock.label, side = 2, line = .5, adj = .5, outer = FALSE, cex = .7) + #mtext("Quarters ahead", side = 1, line = 1.5, adj = .5, outer = FALSE, cex = .7) + + + abline(h=0,lwd = .5) + #abline(v=0) + + } + + options(digits = 4) + + #print the plot + forecasts(ymin= ylim[1], + ymax= ylim[2], + ystep= ystep, + xmin=1, + xmax=nrow(plotData), + xstep=1, + colors=c("firebrick","firebrick","firebrick"), + legcol=1, + keynames=c(""), + hadj=.5, + #note="Note: IRFs are calculated with zero short-run restrictions (cholesky) and represent the response to a one standard deviation \nshock to the Alternative FCI. The blue bands represent the 95 percent confidence interval around the point estimate in red.", + footnote.placement=2.5, + units = ylab, + chart_title = title + #yaxis.shock.label="Percentage Points" + ) + + mtext('Horizon (Quarters)',side = 1, line = .9, adj=.5, cex = .75) + +} diff --git a/R/regimes.R b/R/regimes.R new file mode 100644 index 0000000..3f81454 --- /dev/null +++ b/R/regimes.R @@ -0,0 +1,66 @@ +#------------------------------------------ +# Function to assign regimes via +# unsupervised machine learning +#------------------------------------------ +#' Assign regimes via unsupervised machine learning methods +#' +#' @param data data.frame, matrix, ts, xts, zoo: Endogenous regressors +#' @param regime.n int: number of regimes to estimate (only applies to kmeans) +#' @param engine string: regime assignment technique ('rf' or 'kmeans') +#' +#' @return `data` as a data.frame with a regime column assigning rows to mutually exclusive regimes. If engine = 'rf' is used then regime probabilities will be returned as well. +#' +#' @examples +#' \dontrun{ +#' +#' learn_regimes( +#' data = Data, +#' regime.n = 3, +#' engine = 'kmeans') +#' } +#' +#' @export + +learn_regimes = function( + data, # data.frame, matrix, ts, xts, zoo: Endogenous regressors + regime.n = 2, # int: number of regimes to estimate (only applies to kmeans) + engine = 'rf' # string: regime assignment technique ('rf' or 'kmeans') +){ + + # function warnings + if(!is.matrix(data) & !is.data.frame(data)){ + errorCondition('data must be a matrix or data.frame') + } + if(!is.numeric(regime.n) | regime.n %% 1 != 0 | regime.n <= 0){ + errorCondition('regime.n must be a positive integer') + } + + # cast as data frame if ts, xts, or zoo object + if(is.ts(data) | xts::is.xts(data) | zoo::is.zoo(data)){ + data = data.frame(date = zoo::index(date), data) + } + + # clean data + X = dplyr::select(data, -date) + X.date = data$date + + # assign regimes + if(engine == 'rf'){ + + model = randomForest::randomForest(X) + regime = data.frame(model$votes) + colnames(regime) = paste0('Prob_regime_',c(1:ncol(regime))) + regime$regime = apply(X = regime, MARGIN = 1, which.max) + regime = data.frame(date = X.date, X, regime) + + }else if(engine == 'kmeans'){ + + model = stats::kmeans(X, centers = regime.n) + regime = data.frame(date = X.date, X, model$cluster) + + } + + # return regime assignments + return(regime) + +} diff --git a/R/threshold_var.R b/R/threshold_var.R new file mode 100644 index 0000000..fb2f0ef --- /dev/null +++ b/R/threshold_var.R @@ -0,0 +1,208 @@ + +#------------------------------------------------------------------- +# Function to estimate threshold VAR +# i.e. state-dependent VARs with an exogenous state-variable +#------------------------------------------------------------------- +#' Estimate multi-regime VAR +#' +#' @param data data.frame, matrix, ts, xts, zoo: Endogenous regressors +#' @param regime string: name or regime assignment vector in the design matrix (data) +#' @param p int: lags +#' @param horizon int: forecast horizons +#' @param freq string: frequency of data (day, week, month, quarter, year) +#' +#' @return list of lists, each regime returns its own list with elements `data`, `model`, `forecasts`, `residuals` +#' +#' @examples +#' \dontrun{ +#' threshold_VAR( +#' data = Data, +#' regime = 'regime', +#' p = 1, +#' horizon = 10, +#' freq = 'month') +#' } +#' +#' @export + +# var function +threshold_VAR = function( + data, # data.frame, matrix, ts, xts, zoo: Endogenous regressors + regime, # string: name or regime assignment vector in the design matrix (data) + p = 1, # int: lags + horizon = 10, # int: forecast horizons + freq = 'month' # string: frequency of data (day, week, month, quarter, year) +){ + + # function warnings + if(!is.matrix(data) & !is.data.frame(data)){ + errorCondition('data must be a matrix or data.frame') + } + if(!is.numeric(p) | p %% 1 != 0){ + errorCondition('p must be an integer') + } + if(!is.numeric(horizon) | horizon %% 1 != 0 | horizon <= 0){ + errorCondition('horizon must be a positive integer') + } + if(!freq %in% c('day','week','month','quarter','year')){ + errorCondition("freq must be one of the following strings: 'day','week','month','quarter','year'") + } + if(!regime %in% colnames(data)){ + errorCondition('regime must be the name of a column in data') + } + + # cast as data frame if ts, xts, or zoo object + if(is.ts(data) | xts::is.xts(data) | zoo::is.zoo(data)){ + data = data.frame(date = zoo::index(date), data) + } + + # declare regressors + regressors = colnames(dplyr::select(data, -date, -regime)) + + # create regressors + Y = data.frame(data) %>% + dplyr::select(regressors, date) %>% + n.lag(lags = p) %>% + dplyr::full_join( + dplyr::select(data, regime = regime, date), + by = 'date') + + ### estimate coefficients ---------------------- + + models = Y %>% + # split by regime + dplyr::group_split(regime) %>% + purrr::map(.f = function(Y){ + + # calculate equation by equation + models = as.list(regressors) %>% + purrr::map(.f = function(target){ + + X = Y %>% dplyr::select(dplyr::contains('.l'), target = target) + + # estimate OLS + model = stats::lm(target ~ ., data = X) + + # coefficients + c = broom::tidy(model) %>% dplyr::select(term, coef = estimate) + c$y = target + + se = broom::tidy(model) %>% dplyr::select(term, std.error) + se$y = target + + # return results + return(list(coef = c, se = se)) + }) + + # extract coefficients + coef = + purrr::map(models, .f = function(X){return(X$coef)}) %>% + purrr::reduce(dplyr::bind_rows) %>% + tidyr::pivot_wider(values_from = coef, names_from = term) + + + # extract coefficients + se = + purrr::map(models, .f = function(X){return(X$se)}) %>% + purrr::reduce(dplyr::bind_rows) %>% + tidyr::pivot_wider(values_from = std.error, names_from = term) + + # package for return + model = list(coef = coef, se = se, p = p, freq = freq, horizon = horizon, regime = regime) + + }) + + names(models) = paste0('regime_', unique(Y$regime)) + + ### estimate forecasts ----------------------- + + forecasts = Y %>% + # split by regime + dplyr::group_split(regime) %>% + purrr::map(.f = function(Y){ + + # set appropriate model for the regime + model = models[[paste0('regime_',unique(Y$regime))]] + coef = model$coef + + forecasts = list() + for(i in 1:horizon){ + + # update X + if(i == 1){ + X = Y %>% dplyr::select(dplyr::contains('.l')) + }else{ + X = forecast_prev %>% + n.lag(lags = p) %>% + dplyr::select(dplyr::contains('.l')) + } + + # estimate i-step ahead forecast + forecast = as.matrix(data.frame(1, X)) %*% as.matrix(t(coef[,-1])) + colnames(forecast) = regressors + + # add in dates + forecast = + data.frame( + date = forecast_date( + forecast.date = Y$date, + horizon = i, + freq = freq), + forecast + ) + + # store forecasts + forecasts[[paste0('H_',i)]] = forecast + forecast_prev = forecast + } + + return(forecasts) + + }) + + names(forecasts) = paste0('regime_', unique(Y$regime)) + + # merge forecasts + forecasts = as.list(c(1:horizon)) %>% + purrr::map(.f = function(horizon){ + + r = forecasts %>% + purrr::map(.f = function(regime){ + return(regime[[paste0('H_',horizon)]]) + }) + + r = purrr::reduce(r, dplyr::bind_rows) %>% + dplyr::arrange(date) %>% + dplyr::left_join( + dplyr::select(Y, regime, date), + by = 'date') + + }) + + names(forecasts) = paste0('H_', c(1:horizon)) + + + + ### calculate residuals ----------------------- + + residuals = forecasts %>% + # error by forecast horizon + purrr::map(.f = function(forecast){ + + error = data.frame(forecast) + error[,c(regressors)] = forecast[,c(regressors)] - data.frame(data)[, c(regressors)] + + return(error) + + }) + + ### return output -------------- + return( + list( + model = models, + data = data, + forecasts = forecasts, + residuals = residuals + ) + ) +} diff --git a/R/threshold_var_fevd.R b/R/threshold_var_fevd.R new file mode 100644 index 0000000..1286593 --- /dev/null +++ b/R/threshold_var_fevd.R @@ -0,0 +1,66 @@ +#-------------------------------------------------------- +# Wrapper function to estimate forecast error variance +#-------------------------------------------------------- +#' Estimate multi-regime forecast error variance decomposition +#' +#' @param threshold_var threshold_var output +#' @param horizon int: number of periods +#' +#' @return list, each regime returns its own long-form data.frame +#' +#' @examples +#' \dontrun{ +#' threshold_var_fevd( +#' threshold_var, +#' bootstraps.num = 10) +#' } +#' +#' @export + +# uses the fevd function found in var_fevd.R +threshold_var_fevd = function( + threshold_var, # threshold_VAR output + horizon = 10 # int: number of periods +){ + + # function warnings + if(!is.numeric(horizon) | horizon %% 1 != 0 | horizon <= 0){ + errorCondition('horizon must be a positive integer') + } + + # set data + data = threshold_var$data + regime = threshold_var$model[[1]]$regime + regimes = unlist(unique(dplyr::select(data, regime))) + regressors = colnames(dplyr::select(data, -date, -regime)) + + # estimate impulse responses by regime + results = split(regimes, seq_along(regimes)) %>% + purrr::map(.f = function(regime.val){ + + # set regime specific data + coef = threshold_var$model[[paste0('regime_',regime.val)]]$coef + residuals = threshold_var$residuals[[1]] %>% dplyr::filter(regime == regime.val) + + # error covariance matrix + cov.matrix = var(na.omit(dplyr::select(residuals, -date, -regime))) + + # forecast error variance decomposition + errors = fevd(Phi = coef[,-c(1,2)], Sig = cov.matrix, lag = horizon) + + # reorganize results + response = data.frame(t(errors$OmegaR)) + response$shock = rep(regressors, horizon + 1) + response$horizon = sort(rep(c(0:horizon), length(regressors))) + response = response %>% dplyr::arrange(shock, horizon) + rownames(response) = NULL + + return(response) + + }) + + names(results) = paste0('regime_', regimes) + + return(results) + +} diff --git a/R/threshold_var_irf.R b/R/threshold_var_irf.R new file mode 100644 index 0000000..fdcaf7e --- /dev/null +++ b/R/threshold_var_irf.R @@ -0,0 +1,341 @@ +#------------------------------------------ +# Function to estimate impulse responses +#------------------------------------------ +#' Estimate multi-regime impulse response functions +#' +#' @param threshold_var threshold_VAR output +#' @param horizon int: number of periods +#' @param bootstraps.num int: number of bootstraps +#' @param CI numeric vector: c(lower ci bound, upper ci bound) +#' +#' @return list of lists, each regime returns its own list with elements `irfs`, `ci.lower`, and `ci.upper`; all elements are long-form data.frames +#' +#' @examples +#' \dontrun{ +#' threshold_var_irf( +#' threshold_var, +#' bootstraps.num = 10, +#' CI = c(0.05,0.95)) +#' } +#' +#' @export + +threshold_var_irf = function( + threshold_var, # threshold VAR output + horizon = 10, # int: number of periods + bootstraps.num = 100, # int: number of bootstraps + CI = c(0.1, 0.9) # numeric vector: c(lower ci bound, upper ci bound) +){ + + # function warnings + if(!is.numeric(bootstraps.num) | bootstraps.num %% 1 != 0){ + errorCondition('bootstraps.num must be an integer') + } + if(!is.numeric(CI) | length(CI) != 2 | min(CI) < 0 | max(CI) > 1 | is.na(sum(CI))){ + errorCondition('CI must be a two element numeric vector bound [0,1]') + } + if(!is.numeric(horizon) | horizon %% 1 != 0 | horizon <= 0){ + errorCondition('horizon must be a positive integer') + } + + # set data + residuals = threshold_var$residuals[[1]] + data = threshold_var$data + regime = threshold_var$model[[1]]$regime + + p = threshold_var$model[[1]]$p + freq = threshold_var$model[[1]]$freq + + regressors = colnames(dplyr::select(data, -date, -regime)) + regimes = unlist(unique(dplyr::select(data, regime))) + + p.lower = CI[1] + p.upper = CI[2] + + # estimate impulse responses by regime + results = split(regimes, seq_along(regimes)) %>% + purrr::map(.f = function(regime.val){ + + # set regime specific data + coef = threshold_var$model[[paste0('regime_',regime.val)]]$coef + residuals = residuals %>% + dplyr::filter(regime == regime.val) %>% + dplyr::select(-date, -regime) + + ### calculate impulse responses -------------- + # error covariance matrix + cov.matrix = var(na.omit(residuals)) + # cholesky decomposition + cholesky.matrix = t(chol(cov.matrix)) + + # extract responses + irfs = regressors %>% + purrr::map(.f = function(shock){ + + # loop overhead + irf = matrix(ncol = length(regressors), nrow = (horizon+1)) + colnames(irf) = regressors + irf[1,] = cholesky.matrix[,shock] + + # responses per horizon + for(j in 1:(horizon)){ + + # initial impact + if(j == 1){ + + AA = c(cholesky.matrix[,shock], rep(0, length(regressors)*(p)-length(regressors))) + response = AA %*% as.matrix(t(coef[,-c(1,2)])) + + # recursively forecasted impact + }else{ + + require.length = length(AA) + current = as.vector(t(irf[j:j-p+1,])) + A = c(current, rep(0, require.length - length(current))) + response = A %*% as.matrix(t(coef[,-c(1,2)])) + + } + + # store + irf[j+1,] = response + + } + + irf = data.frame(shock = shock, horizon = c(0:horizon), irf) + + # return irf + return(irf) + + }) + + irfs = purrr::reduce(irfs, dplyr::bind_rows) + + ### bootstrap irf standard errors -------------- + # see Lutkepohl (2005) + + # 1. create bootstrap time series + bagged.series = as.list(1:bootstraps.num) %>% + purrr::map(.f = function(count){ + + # draw bootstrapped residuals + U = na.omit(residuals)[sample( + c(1:nrow(na.omit(residuals))), + size = nrow(data), + replace = TRUE), + ] + U = U %>% + dplyr::mutate_all(function(X){return(X-mean(X, na.rm = T))}) + + # create lags + X = data %>% + dplyr::select(-regime) %>% + n.lag(lags = p) %>% + dplyr::select(dplyr::contains('.l')) + + # estimate time series + Y = as.matrix(data.frame(1, X)) %*% as.matrix(t(coef[,-1])) + Y = Y + U + colnames(Y) = regressors + Y = data.frame(Y, date = data$date) + + # filter for regimes + Y = Y %>% + dplyr::full_join( + dplyr::select(data, regime = regime, date), + by = 'date') %>% + dplyr::filter(regime == regime.val) %>% + dplyr::select(-regime) %>% + na.omit() + + # return synthetic observations + return(Y) + + }) + + # 2. create bootstrapped residuals + bagged.irf = bagged.series %>% + purrr::map(.f = function(synth){ + + # re-estimate VAR with bagged series + var.boot = + suppressMessages( + VAR(data = synth, + p = p, + horizon = 1, + freq = freq) + ) + + residuals = var.boot$residuals[[1]] %>% dplyr::select(-date) + coef = var.boot$model$coef + + # error covariance matrix + cov.matrix = var(na.omit(residuals)) + # cholesky decomposition + cholesky.matrix = t(chol(cov.matrix)) + colnames(cholesky.matrix) = regressors + + # extract responses + irfs = regressors %>% + purrr::map(.f = function(shock){ + + # loop overhead + irf = matrix(ncol = length(regressors), nrow = (horizon+1)) + colnames(irf) = regressors + irf[1,] = cholesky.matrix[,shock] + + # responses per horizon + for(j in 1:(horizon)){ + + # initial impact + if(j == 1){ + + AA = c(cholesky.matrix[,shock], rep(0, length(regressors)*(p)-length(regressors))) + response = AA %*% as.matrix(t(coef[,-c(1,2)])) + + # recursively forecasted impact + }else{ + + require.length = length(AA) + current = as.vector(t(irf[j:j-p+1,])) + A = c(current, rep(0, require.length - length(current))) + response = A %*% as.matrix(t(coef[,-c(1,2)])) + + } + + # store + irf[j+1,] = response + + } + + irf = data.frame(shock = shock, horizon = c(0:horizon), irf) + + # return irf + return(irf) + + }) + + irfs = purrr::reduce(irfs, dplyr::bind_rows) + + return(irfs) + + }) + + # 3. calculate confidence intervals + ci.lower = bagged.irf %>% + purrr::reduce(dplyr::bind_rows) %>% + dplyr::group_by(shock, horizon) %>% + dplyr::summarise_all(quantile, p.lower, na.rm = T) %>% + dplyr::arrange(shock, horizon) + + ci.upper = bagged.irf %>% + purrr::reduce(dplyr::bind_rows) %>% + dplyr::group_by(shock, horizon) %>% + dplyr::summarise_all(quantile, p.upper, na.rm = T) %>% + dplyr::arrange(shock, horizon) + + ci.med = bagged.irf %>% + purrr::reduce(dplyr::bind_rows) %>% + dplyr::group_by(shock, horizon) %>% + dplyr::summarise_all(quantile, 0.5, na.rm = T) %>% + dplyr::arrange(shock, horizon) + + irfs = irfs %>% + dplyr::arrange(shock, horizon) + + ci.adjust = irfs[,-c(1,2)] - ci.med[,-c(1,2)] + ci.lower[,-c(1,2)] = ci.lower[,-c(1,2)] + ci.adjust + ci.upper[,-c(1,2)] = ci.upper[,-c(1,2)] + ci.adjust + + irfs = list(irfs = irfs, ci.upper = ci.upper, ci.lower = ci.lower) + + return(irfs) + + }) + + names(results) = paste0('regime_', regimes) + + ### return output -------------- + return(results) + +} + +#------------------------------------------ +# Function to plot IRfs +#------------------------------------------ +### Function to plot individual irf plot +individual_irf_plot = function( + irfs, # var_irf object + shock.var, # string: name of variable to treat as the shock + response.var, # string: name of variable to treat as the response + title, # string: title of the chart + ylab # string: y-axis label +){ + + # set data + irf = irfs$irfs + irf.lower = irfs$ci.lower + irf.upper = irfs$ci.upper + + # filter for one shock + irf = irf %>% filter(shock == shock.var) + irf.lower = irf.lower %>% filter(shock == shock.var) + irf.upper = irf.upper %>% filter(shock == shock.var) + + # filter for one response + irf = irf %>% select(point = response.var, horizon) + irf.lower = irf.lower %>% ungroup() %>% select(lower = response.var) + irf.upper = irf.upper %>% ungroup() %>% select(upper = response.var) + + plotdata = cbind(irf.lower, irf, irf.upper) + + # plot GDP + response.gdp <- plotdata %>% + ggplot(aes(x=horizon, y=point, ymin=lower, ymax=upper)) + + geom_hline(yintercept = 0, color="red") + + geom_ribbon(fill="grey", alpha=0.2) + + geom_line() + + theme_light() + + ggtitle(title)+ + ylab(ylab)+ + xlab("") + + theme(plot.title = element_text(size = 11, hjust=0.5), + axis.title.y = element_text(size=11)) + + response.gdp + +} + +### function to plot all irfs +irf_plot = function( + irfs, # var_irf object + shocks, # string vector: shocks to plot + responses # string vector: responses to plot +){ + + # set shocks and responses + shocks = unique(irfs$irfs$shock) + responses = unique(irfs$irfs$shock) + + # generate plots + plot.names = expand_grid(shock = shocks, response = responses) + plots = split(plot.names, seq(nrow(plot.names))) %>% + purrr::map(.f = function(x){ + + chart = + individual_irf_plot( + irfs, + shock.var = x$shock, + response.var = x$response, + title = paste0(x$response, ' response to ', x$shock, ' shock'), + ylab = '') + + return(chart) + + }) + + # create plot + n <- length(plots) + nCol <- floor(sqrt(n)) + do.call(gridExtra::grid.arrange, c(plots, ncol=nCol)) + +} diff --git a/R/var.R b/R/var.R new file mode 100644 index 0000000..8e695ef --- /dev/null +++ b/R/var.R @@ -0,0 +1,152 @@ +#------------------------------------------ +# Function to estimate VAR +#------------------------------------------ +#' Estimate single-regime VAR +#' +#' @param data data.frame, matrix, ts, xts, zoo: Endogenous regressors +#' @param p int: lags +#' @param horizon int: forecast horizons +#' @param freq string: frequency of data (day, week, month, quarter, year) +#' +#' @return list object with elements `data`, `model`, `forecasts`, `residuals` +#' +#' @examples +#' \dontrun{ +#' VAR( +#' data = Data, +#' p = 1, +#' horizon = 10, +#' freq = 'month') +#' } +#' +#' @export + +# var function +VAR = function( + data, # data.frame, matrix, ts, xts, zoo: Endogenous regressors + p = 1, # int: lags + horizon = 10, # int: forecast horizons + freq = 'month' # string: frequency of data (day, week, month, quarter, year) +){ + + # function warnings + if(!is.matrix(data) & !is.data.frame(data)){ + errorCondition('data must be a matrix or data.frame') + } + if(!is.numeric(p) | p %% 1 != 0){ + errorCondition('p must be an integer') + } + if(!is.numeric(horizon) | horizon %% 1 != 0 | horizon <= 0){ + errorCondition('horizon must be a positive integer') + } + if(!freq %in% c('day','week','month','quarter','year')){ + errorCondition("freq must be one of the following strings: 'day','week','month','quarter','year'") + } + + # cast as data frame if ts, xts, or zoo object + if(is.ts(data) | xts::is.xts(data) | zoo::is.zoo(data)){ + data = data.frame(date = zoo::index(date), data) + } + + # declare regressors + regressors = colnames(dplyr::select(data, -date)) + + # create regressors + Y = data.frame(data) %>% + n.lag(lags = p) + + # remove date + Y = Y %>% dplyr::select(-date) + + + ### estimate coefficients ---------------------- + models = as.list(regressors) %>% + purrr::map(.f = function(target){ + + X = Y %>% dplyr::select(dplyr::contains('.l'), target = target) + + # estimate OLS + model = stats::lm(target ~ ., data = X) + + # coefficients + c = broom::tidy(model) %>% dplyr::select(term, coef = estimate) + c$y = target + + se = broom::tidy(model) %>% dplyr::select(term, std.error) + se$y = target + + # return results + return(list(coef = c, se = se)) + }) + + # extract coefficients + coef = + purrr::map(models, .f = function(X){return(X$coef)}) %>% + purrr::reduce(dplyr::bind_rows) %>% + tidyr::pivot_wider(values_from = coef, names_from = term) + + + # extract coefficients + se = + purrr::map(models, .f = function(X){return(X$se)}) %>% + purrr::reduce(dplyr::bind_rows) %>% + tidyr::pivot_wider(values_from = std.error, names_from = term) + + # package for return + model = list(coef = coef, se = se, p = p, freq = freq, horizon = horizon) + + ### estimate forecasts ----------------------- + forecasts = list() + for(i in 1:horizon){ + + # update X + if(i == 1){ + X = Y %>% dplyr::select(dplyr::contains('.l')) + }else{ + X = forecast_prev %>% + n.lag(lags = p) %>% + dplyr::select(dplyr::contains('.l')) + } + + # estimate i-step ahead forecast + forecast = as.matrix(data.frame(1, X)) %*% as.matrix(t(coef[,-1])) + colnames(forecast) = regressors + + # add in dates + forecast = + data.frame( + date = forecast_date( + forecast.date = data$date, + horizon = i, + freq = freq), + forecast + ) + + # store forecasts + forecasts[[paste0('H_',i)]] = forecast + forecast_prev = forecast + } + + ### calculate residuals ----------------------- + residuals = forecasts %>% + # error by forecast horizon + purrr::map(.f = function(forecast){ + + error = data.frame(forecast) + error[,c(regressors)] = forecast[,c(regressors)] - data.frame(data)[, c(regressors)] + + return(error) + + }) + + + ### return output -------------- + return( + list( + model = model, + data = data, + forecasts = forecasts, + residuals = residuals + ) + ) +} diff --git a/R/var_fevd.R b/R/var_fevd.R new file mode 100644 index 0000000..676114a --- /dev/null +++ b/R/var_fevd.R @@ -0,0 +1,142 @@ +#--------------------------------------------- +# Estimate forecast error variance +# source code adapted from the MTS package +#--------------------------------------------- +fevd = function(Phi, Sig, lag = 4) +{ + if (length(Phi) > 0) { + if (!is.matrix(Phi)) + Phi = as.matrix(Phi) + } + + if (!is.matrix(Sig)) + Sig = as.matrix(Sig) + if (lag < 1) + lag = 1 + p = 0 + if (length(Phi) > 0) { + k = nrow(Phi) + m = ncol(Phi) + p = floor(m/k) + } + q = 0 + + Si = diag(rep(1, k)) + + m = (lag + 1) * k + m1 = (q + 1) * k + if (m > m1) { + Si = cbind(Si, matrix(0, k, (m - m1))) + } + if (p > 0) { + for (i in 1:lag) { + if (i <= p) { + idx = (i - 1) * k + tmp = Phi[, (idx + 1):(idx + k)] + } + else { + tmp = matrix(0, k, k) + } + jj = i - 1 + jp = min(jj, p) + if (jp > 0) { + for (j in 1:jp) { + jdx = (j - 1) * k + idx = (i - j) * k + w1 = Phi[, (jdx + 1):(jdx + k)] + w2 = Si[, (idx + 1):(idx + k)] + tmp = tmp + w1 %*% w2 + } + } + kdx = i * k + Si[, (kdx + 1):(kdx + k)] = tmp + } + } + orSi = NULL + m1 = chol(Sig) + P = t(m1) + orSi = P + for (i in 1:lag) { + idx = i * k + w1 = Si[, (idx + 1):(idx + k)] + w2 = w1 %*% P + orSi = cbind(orSi, w2) + } + orSi2 = orSi^2 + Ome = orSi2[, 1:k] + wk = Ome + for (i in 1:lag) { + idx = i * k + wk = wk + orSi2[, (idx + 1):(idx + k)] + Ome = cbind(Ome, wk) + } + FeV = NULL + OmeRa = Ome[, 1:k] + FeV = cbind(FeV, apply(OmeRa, 1, sum)) + OmeRa = OmeRa/FeV[, 1] + for (i in 1:lag) { + idx = i * k + wk = Ome[, (idx + 1):(idx + k)] + FeV = cbind(FeV, apply(wk, 1, sum)) + OmeRa = cbind(OmeRa, wk/FeV[, (i + 1)]) + } + for (i in 1:(lag + 1)) { + idx = (i - 1) * k + Ratio = OmeRa[, (idx + 1):(idx + k)] + } + FEVdec <- list(irf = Si, orthirf = orSi, Omega = Ome, OmegaR = OmeRa) + + return(FEVdec) +} + +#-------------------------------------------------------- +# Wrapper function to estimate forecast error variance +#-------------------------------------------------------- +#' Estimate single-regime forecast error variance decomposition +#' +#' @param var VAR output +#' @param horizon int: number of periods +#' +#' @return long-form data.frame +#' +#' @examples +#' \dontrun{ +#' var_fevd( +#' var, +#' bootstraps.num = 10) +#' } +#' +#' @export + +var_fevd = function( + var, # VAR output + horizon = 10 # int: number of periods +){ + + # function warnings + if(!is.numeric(horizon) | horizon %% 1 != 0 | horizon <= 0){ + errorCondition('horizon must be a positive integer') + } + + # set data + coef = var$model$coef + residuals = var$residuals[[1]] + data = var$data + regressors = colnames(dplyr::select(data, -date)) + + # error covariance matrix + cov.matrix = var(na.omit(dplyr::select(residuals, -date))) + + # forecast error variance decomposition + errors = fevd(Phi = coef[,-c(1,2)], Sig = cov.matrix, lag = horizon) + + # reorganize results + response = data.frame(t(errors$OmegaR)) + response$shock = rep(regressors, horizon + 1) + response$horizon = sort(rep(c(0:horizon), length(regressors))) + response = response %>% dplyr::arrange(shock, horizon) + rownames(response) = NULL + + return(response) + +} diff --git a/R/var_irf.R b/R/var_irf.R new file mode 100644 index 0000000..2caf426 --- /dev/null +++ b/R/var_irf.R @@ -0,0 +1,340 @@ +#------------------------------------------ +# Function to estimate impulse responses +#------------------------------------------ +#' Estimate single-regime impulse response functions +#' +#' @param var VAR output +#' @param horizon int: number of periods +#' @param bootstraps.num int: number of bootstraps +#' @param CI numeric vector: c(lower ci bound, upper ci bound) +#' +#' @return list object with elements `irfs`, `ci.lower`, and `ci.upper`; all elements are long-form data.frames +#' +#' @examples +#' \dontrun{ +#' var_irf( +#' var, +#' bootstraps.num = 10, +#' CI = c(0.05,0.95)) +#' } +#' +#' @export + +var_irf = function( + var, # VAR output + horizon = 10, # int: number of periods + bootstraps.num = 100, # int: number of bootstraps + CI = c(0.1, 0.9) # numeric vector: c(lower ci bound, upper ci bound) +){ + + # function warnings + if(!is.numeric(bootstraps.num) | bootstraps.num %% 1 != 0){ + errorCondition('bootstraps.num must be an integer') + } + if(!is.numeric(CI) | length(CI) != 2 | min(CI) < 0 | max(CI) > 1 | is.na(sum(CI))){ + errorCondition('CI must be a two element numeric vector bound [0,1]') + } + if(!is.numeric(horizon) | horizon %% 1 != 0 | horizon <= 0){ + errorCondition('horizon must be a positive integer') + } + + # set data + coef = var$model$coef + residuals = var$residuals + data = var$data + + # set variables + p = var$model$p + freq = var$model$freq + horizon = length(residuals) + regressors = colnames(dplyr::select(data, -date)) + p.lower = CI[1] + p.upper = CI[2] + + ### calculate impulse responses -------------- + # error covariance matrix + cov.matrix = var(na.omit(dplyr::select(residuals[[1]], -date))) + # cholesky decomposition + cholesky.matrix = t(chol(cov.matrix)) + + # extract responses + irfs = regressors %>% + purrr::map(.f = function(shock){ + + # loop overhead + irf = matrix(ncol = length(regressors), nrow = (horizon+1)) + colnames(irf) = regressors + irf[1,] = cholesky.matrix[,shock] + + # responses per horizon + for(j in 1:(horizon)){ + + # initial impact + if(j == 1){ + + AA = c(cholesky.matrix[,shock], rep(0, length(regressors)*(p)-length(regressors))) + response = AA %*% as.matrix(t(coef[,-c(1,2)])) + + # recursively forecasted impact + }else{ + + require.length = length(AA) + current = as.vector(t(irf[j:j-p+1,])) + A = c(current, rep(0, require.length - length(current))) + response = A %*% as.matrix(t(coef[,-c(1,2)])) + + } + + # store + irf[j+1,] = response + + } + + irf = data.frame(shock = shock, horizon = c(0:horizon), irf) + + # return irf + return(irf) + + }) + + irfs = purrr::reduce(irfs, dplyr::bind_rows) + + ### bootstrap irf standard errors -------------- + # see Lutkepohl (2005) + + # 1. create bootstrap time series + bagged.series = as.list(1:bootstraps.num) %>% + purrr::map(.f = function(count){ + + # draw bootstrapped residuals + U = residuals[[1]][sample(c(1:nrow(residuals[[1]])), + size = nrow(residuals[[1]]), + replace = TRUE),] + U = U %>% + dplyr::select(-date) %>% + dplyr::mutate_all(function(X){return(X-mean(X, na.rm = T))}) + + # create lags + X = data %>% + n.lag(lags = p) %>% + dplyr::select(dplyr::contains('.l')) + + # estimate time series + Y = as.matrix(data.frame(1, X)) %*% as.matrix(t(coef[,-1])) + Y = Y + U + colnames(Y) = regressors + Y = data.frame(Y, date = data$date) + + # return synthetic observations + return(Y) + + }) + + # 2. create bootstrapped residuals + bagged.irf = bagged.series %>% + purrr::map(.f = function(synth){ + + # re-estimate VAR with bagged series + var.boot = + suppressMessages( + VAR(data = synth, + p = p, + horizon = 1, + freq = freq) + ) + + residuals = var.boot$residuals + coef = var.boot$model$coef + + # error covariance matrix + cov.matrix = var(na.omit(dplyr::select(residuals[[1]], -date))) + # cholesky decomposition + cholesky.matrix = t(chol(cov.matrix)) + colnames(cholesky.matrix) = regressors + + # extract responses + irfs = regressors %>% + purrr::map(.f = function(shock){ + + # loop overhead + irf = matrix(ncol = length(regressors), nrow = (horizon+1)) + colnames(irf) = regressors + irf[1,] = cholesky.matrix[,shock] + + # responses per horizon + for(j in 1:(horizon)){ + + # initial impact + if(j == 1){ + + AA = c(cholesky.matrix[,shock], rep(0, length(regressors)*(p)-length(regressors))) + response = AA %*% as.matrix(t(coef[,-c(1,2)])) + + # recursively forecasted impact + }else{ + + require.length = length(AA) + current = as.vector(t(irf[j:j-p+1,])) + A = c(current, rep(0, require.length - length(current))) + response = A %*% as.matrix(t(coef[,-c(1,2)])) + + } + + # store + irf[j+1,] = response + + } + + irf = data.frame(shock = shock, horizon = c(0:horizon), irf) + + # return irf + return(irf) + + }) + + irfs = purrr::reduce(irfs, dplyr::bind_rows) + + return(irfs) + + }) + + # 3. calculate confidence intervals + ci.lower = bagged.irf %>% + purrr::reduce(dplyr::bind_rows) %>% + dplyr::group_by(shock, horizon) %>% + dplyr::summarise_all(quantile, p.lower, na.rm = T) %>% + dplyr::arrange(shock, horizon) + + ci.upper = bagged.irf %>% + purrr::reduce(dplyr::bind_rows) %>% + dplyr::group_by(shock, horizon) %>% + dplyr::summarise_all(quantile, p.upper, na.rm = T) %>% + dplyr::arrange(shock, horizon) + + ci.med = bagged.irf %>% + purrr::reduce(dplyr::bind_rows) %>% + dplyr::group_by(shock, horizon) %>% + dplyr::summarise_all(quantile, 0.5, na.rm = T) %>% + dplyr::arrange(shock, horizon) + + irfs = irfs %>% + dplyr::arrange(shock, horizon) + + ci.adjust = irfs[,-c(1,2)] - ci.med[,-c(1,2)] + ci.lower[,-c(1,2)] = ci.lower[,-c(1,2)] + ci.adjust + ci.upper[,-c(1,2)] = ci.upper[,-c(1,2)] + ci.adjust + + irfs = list(irfs = irfs, ci.upper = ci.upper, ci.lower = ci.lower) + + ### return output -------------- + return(irfs) +} + + + + +#------------------------------------------ +# Function to plot IRfs +#------------------------------------------ +### Function to plot individual irf plot + +#' Plot an individual IRF +#' +#' @param irfs var_irf object +#' @param shock.var string: name of variable to treat as the shock +#' @param response.var string: name of variable to treat as the response +#' @param title string: title of the chart +#' @param ylab string: y-axis label +#' +#' @return ggplot2 graph +#' +#' @export +individual_var_irf_plot = function( + irfs, # var_irf object + shock.var, # string: name of variable to treat as the shock + response.var, # string: name of variable to treat as the response + title, # string: title of the chart + ylab # string: y-axis label +){ + + # set data + irf = irfs$irfs + irf.lower = irfs$ci.lower + irf.upper = irfs$ci.upper + + # filter for one shock + irf = irf %>% filter(shock == shock.var) + irf.lower = irf.lower %>% filter(shock == shock.var) + irf.upper = irf.upper %>% filter(shock == shock.var) + + # filter for one response + irf = irf %>% select(point = response.var, horizon) + irf.lower = irf.lower %>% ungroup() %>% select(lower = response.var) + irf.upper = irf.upper %>% ungroup() %>% select(upper = response.var) + + plotdata = cbind(irf.lower, irf, irf.upper) + + # plot GDP + response.gdp <- plotdata %>% + ggplot(aes(x=horizon, y=point, ymin=lower, ymax=upper)) + + geom_hline(yintercept = 0, color="red") + + geom_ribbon(fill="grey", alpha=0.2) + + geom_line() + + theme_light() + + ggtitle(title)+ + ylab(ylab)+ + xlab("") + + theme(plot.title = element_text(size = 11, hjust=0.5), + axis.title.y = element_text(size=11)) + + response.gdp + +} + +### function to plot all irfs +#' Plot all IRFs +#' +#' @param irfs var_irf object +#' @param shocks string vector: shocks to plot +#' @param responses string vector: responses to plot +#' +#' @return grid of ggplot2 graphs +#' +#' @export + +var_irf_plot = function( + irfs, # var_irf object + shocks, # string vector: shocks to plot + responses # string vector: responses to plot +){ + + # function variables + shock = repsonse = NA + + # set shocks and responses + shocks = unique(irfs$irfs$shock) + responses = unique(irfs$irfs$shock) + + # generate plots + plot.names = tidyr::expand_grid(shock = shocks, response = responses) + plots = split(plot.names, seq(nrow(plot.names))) %>% + purrr::map(.f = function(x){ + + chart = + individual_irf_plot( + irfs, + shock.var = x$shock, + response.var = x$response, + title = paste0(x$response, ' response to ', x$shock, ' shock'), + ylab = '') + + return(chart) + + }) + + # create plot + n <- length(plots) + nCol <- floor(sqrt(n)) + do.call(gridExtra::grid.arrange, c(plots, ncol=nCol)) + +} diff --git a/README.md b/README.md new file mode 100644 index 0000000..0e19e2f --- /dev/null +++ b/README.md @@ -0,0 +1,162 @@ +# sovereign: State-Dependent Empirical Analysis + + +[![License: GPL v3](https://img.shields.io/badge/License-GPL%20v3-blue.svg)](http://www.gnu.org/licenses/gpl-3.0) +[![Lifecycle: experimental](https://img.shields.io/badge/lifecycle-experimental-orange.svg)](https://www.tidyverse.org/lifecycle/#experimental) +[![R-CMD-check](https://github.com/r-lib/usethis/workflows/R-CMD-check/badge.svg)](https://github.com/r-lib/usethis/actions) + + +The `sovereign` package introduces a set of tools for state-dependant empirical analysis through both VAR- and local projection-based state-dependant forecasts, impulse response functions, and forecast error variance decomposition. + +The `sovereign` package remains under active development. As a result, **the API is not to be considered stable**, and future updates will most likely deprecate and break current functions. + +## Available Tools + +Unsupervised Regime Assignment +1. random forest +2. k-means clustering + +Local Projections +1. impulse responses & charting + +Vector Auto-Regression (VAR) +1. recursive forecasting +2. forecast error variance decomposition +3. impulse responses with bootstrapped confidence intervals and charting + +---- + +## Basic Workflow + # load packages + library(sovereign) # analysis + library(tidyverse) # general cleaning + library(lubridate) # date functions + + #------------------------------------------- + # create data + #------------------------------------------- + # pull and prepare data from FRED + quantmod::getSymbols.FRED( + c('UNRATE','INDPRO','GS10'), + env = globalenv()) + + Data = cbind(UNRATE, INDPRO, GS10) + + Data = data.frame(Data, date = zoo::index(Data)) %>% + dplyr::filter(lubridate::year(date) >= 1990) %>% + na.omit() + + # create a regime explicitly + Data.threshold = Data %>% + mutate(mp = if_else(GS10 > median(GS10), 1, 0)) + + #------------------------------------------ + # learn regimes + #------------------------------------------ + # assign regimes based on unsurpervised kmeans + # (will not be used further) + regimes = + learn_regimes( + data = Data, + regime.n = 3, + engine = 'kmeans') + + #------------------------------------------ + # single-regime var + #------------------------------------------ + # estimate VAR + var = + VAR( + data = Data, + p = 1, + horizon = 10, + freq = 'month') + + # estimate IRF + irf = + var_irf_plot( + var, + bootstraps.num = 10, + CI = c(0.05,0.95)) + + # plot IRF + irf_plot(irf) + + # estimate forecast error variance decomposition + fevd = + var_fevd( + var, + horizon = 10) + + #------------------------------------------- + # multi-regime var + #------------------------------------------- + # estimate multi-regime VAR + tvar = + threshold_VAR( + data = Data.threshold, + regime = 'mp', + p = 1, + horizon = 1, + freq = 'month') + + # estimate IRF + tvar.irf = + threshold_var_irf( + tvar, + horizon = 10, + bootstraps.num = 10, + CI = c(0.05,0.95)) + + # plot IRF + # regime 1: low interest rates + irf_plot(tvar.irf[[1]]) + # regime 2: high interest rates + irf_plot(tvar.irf[[2]]) + + # estimate forecast error variance decomposition + tvar.fevd = + threshold_var_fevd( + tvar, + horizon = 10) + + #------------------------------------------- + # local projection IRFs + #------------------------------------------- + # estimate single-regime IRF + lp.irf = + lp_irf( + data = Data, + shock = 'INDPRO', + target = 'GS10', + horizons = 20, + lags = 2) + + # estimate multi-regime IRF + tlp.irf = + threshold_lp_irf( + data = dplyr::select(Data.threshold, -reg), + thresholdVar = dplyr::select(Data, reg, date), + shock = 'AA', + target = 'BB', + horizons = 20, + lags = 2) + + + +--- +## Known problems / wishlist +Code +1. local projection forecasting +2. tvar forecasting is restricted to one period ahead +3. add confidence intervals to forecast error variance decomposition +4. clean local projection functions +5. plot tests + +Package +1. ~~documentation~~ +2. ~~tests~~ +3. ~~simple vignette~~ +4. ~~badges~~ +4. ~~github~~ +5. add plotting to example diff --git a/codecov.yml b/codecov.yml new file mode 100644 index 0000000..2ae0ff7 --- /dev/null +++ b/codecov.yml @@ -0,0 +1,19 @@ +comment: false +language: R +sudo: false +cache: packages +after_success: +- Rscript -e 'covr::codecov()' + +coverage: + status: + project: + default: + target: auto + threshold: 1% + informational: true + patch: + default: + target: auto + threshold: 1% + informational: true diff --git a/development_tests.R b/development_tests.R new file mode 100644 index 0000000..bba23e3 --- /dev/null +++ b/development_tests.R @@ -0,0 +1,108 @@ + +############################################################################################################### +# Function test +############################################################################################################### +setwd('/scratch/m1tjp01/Tyler/sovereign/') + +# ThresholdVAR files +source('./R/helper.R') +source('./R/var_fevd.R') +source('./R/var_irf.R') +source('./R/var.R') +source('./R/threshold_var_fevd.R') +source('./R/threshold_var_irf.R') +source('./R/threshold_var.R') + +# libraries +library(tis) +library(stfm.helper, lib.loc="/stfm/shared1/R") +library(tidyverse) # general cleaning +library(lubridate) # date functions + +#------------------------------------------- +# create data +#------------------------------------------ +# load in unemployment and gdp +usTickers = c('gdp_xcw_09.q','ruc.q') +econActivity = getfame_dt(usTickers,"us") %>% + rename(rgdp = gdp_xcw_09.q, unemp = ruc.q) %>% + mutate(unemp = unemp - lag(unemp)) +day(econActivity$date) <- 1 + +# load effective fed funds data +ff = as.data.frame(getfame_dt("rifspff_n.b","us")) %>% + group_by(year = year(date), quarter = quarter(date)) %>% + summarise(fedfunds = mean(rifspff_n.b, na.rm = T)) %>% + mutate(date = lubridate::ymd(paste0(year,'-',(quarter*3),'-01'))) %>% + ungroup() %>% select(date,fedfunds) + +Data = full_join(econActivity, ff, by = 'date') %>% na.omit() + +#------------------------------------------- +# test single-regime var +#------------------------------------------ +# estimate VAR +test.var = + VAR( + data = Data, + p = 1, + horizon = 10, + freq = 'quarter') + +# estimate IRF +test.irf = + var_irf(test.var, + bootstraps.num = 10, + CI = c(0.05,0.95)) + +# plot IRF +irf_plot(test.irf) + +# forecast error variance decomposition +test.fevd = + var_fevd( + var = test.var, + horizon = 10) + +#------------------------------------------- +# baseline testing against vars package +#------------------------------------------ +# test.var coefficients align with the baseline.var coef! +baseline.var = + vars::VAR(y = select(Data, -date), p = 1, type = 'const') + +# test.irf coefficients align with the baseline.irf coef! +baseline.irf = + vars::irf(baseline.var) + +# exact match +baseline.fevd = + vars::fevd(baseline.var) + +#------------------------------------------- +# test multi-regime var +#------------------------------------------ + +Data.threshold = Data %>% + mutate(mp = if_else(fedfunds > median(fedfunds), 1, 0)) + +test.tvar = + threshold_VAR( + data = Data.threshold, + regime = 'mp', + p = 1, + horizon = 1, + freq = 'quarter' + ) + +test.tvar.irf = + threshold_var_irf( + test.tvar, + horizon = 10, + bootstraps.num = 10, + CI = c(0.05,0.95)) + +test.tvar.fevd = + threshold_var_fevd( + test.tvar, + horizon = 10) diff --git a/man/VAR.Rd b/man/VAR.Rd new file mode 100644 index 0000000..2a84252 --- /dev/null +++ b/man/VAR.Rd @@ -0,0 +1,33 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/var.R +\name{VAR} +\alias{VAR} +\title{Estimate single-regime VAR} +\usage{ +VAR(data, p = 1, horizon = 10, freq = "month") +} +\arguments{ +\item{data}{data.frame, matrix, ts, xts, zoo: Endogenous regressors} + +\item{p}{int: lags} + +\item{horizon}{int: forecast horizons} + +\item{freq}{string: frequency of data (day, week, month, quarter, year)} +} +\value{ +list object with elements \code{data}, \code{model}, \code{forecasts}, \code{residuals} +} +\description{ +Estimate single-regime VAR +} +\examples{ +\dontrun{ +VAR( + data = Data, + p = 1, + horizon = 10, + freq = 'month') +} + +} diff --git a/man/individual_var_irf_plot.Rd b/man/individual_var_irf_plot.Rd new file mode 100644 index 0000000..1e971d3 --- /dev/null +++ b/man/individual_var_irf_plot.Rd @@ -0,0 +1,25 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/var_irf.R +\name{individual_var_irf_plot} +\alias{individual_var_irf_plot} +\title{Plot an individual IRF} +\usage{ +individual_var_irf_plot(irfs, shock.var, response.var, title, ylab) +} +\arguments{ +\item{irfs}{var_irf object} + +\item{shock.var}{string: name of variable to treat as the shock} + +\item{response.var}{string: name of variable to treat as the response} + +\item{title}{string: title of the chart} + +\item{ylab}{string: y-axis label} +} +\value{ +ggplot2 graph +} +\description{ +Plot an individual IRF +} diff --git a/man/learn_regimes.Rd b/man/learn_regimes.Rd new file mode 100644 index 0000000..3aa883e --- /dev/null +++ b/man/learn_regimes.Rd @@ -0,0 +1,31 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/regimes.R +\name{learn_regimes} +\alias{learn_regimes} +\title{Assign regimes via unsupervised machine learning methods} +\usage{ +learn_regimes(data, regime.n = 2, engine = "rf") +} +\arguments{ +\item{data}{data.frame, matrix, ts, xts, zoo: Endogenous regressors} + +\item{regime.n}{int: number of regimes to estimate (only applies to kmeans)} + +\item{engine}{string: regime assignment technique ('rf' or 'kmeans')} +} +\value{ +\code{data} as a data.frame with a regime column assigning rows to mutually exclusive regimes. If engine = 'rf' is used then regime probabilities will be returned as well. +} +\description{ +Assign regimes via unsupervised machine learning methods +} +\examples{ +\dontrun{ + + learn_regimes( + data = Data, + regime.n = 3, + engine = 'kmeans') +} + +} diff --git a/man/lp_irf.Rd b/man/lp_irf.Rd new file mode 100644 index 0000000..44b65c2 --- /dev/null +++ b/man/lp_irf.Rd @@ -0,0 +1,36 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/local_projections.R +\name{lp_irf} +\alias{lp_irf} +\title{Estimate single-regime local projection IRFs} +\usage{ +lp_irf(data, shock, target, horizons, lags) +} +\arguments{ +\item{data}{data.frame, matrix, ts, xts, zoo: Endogenous regressors} + +\item{shock}{string: variable to shock} + +\item{target}{string: variable betas to collect} + +\item{horizons}{int: horizons to forecast out to} + +\item{lags}{int: lags to include in regressions} +} +\value{ +data.frame +} +\description{ +Estimate single-regime local projection IRFs +} +\examples{ +\dontrun{ +lp_irf( + data = Data, + shock = 'x', + target = 'y', + horizons = 20, + lags = 2) +} + +} diff --git a/man/lp_irf_chart.Rd b/man/lp_irf_chart.Rd new file mode 100644 index 0000000..51b53f8 --- /dev/null +++ b/man/lp_irf_chart.Rd @@ -0,0 +1,25 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/lp_irf_chart.R +\name{lp_irf_chart} +\alias{lp_irf_chart} +\title{Chart local projection IRFs} +\usage{ +lp_irf_chart(plotData, title, ylim, ystep, ylab) +} +\arguments{ +\item{plotData}{lp_irf output} + +\item{title}{string: chart title} + +\item{ylim}{int: y-axis limits, c(lower limit, upper limit)} + +\item{ystep}{int: step size inbetween y-axis tick marks} + +\item{ylab}{string: y-axis label} +} +\value{ +plot +} +\description{ +Chart local projection IRFs +} diff --git a/man/pipe.Rd b/man/pipe.Rd new file mode 100644 index 0000000..62311aa --- /dev/null +++ b/man/pipe.Rd @@ -0,0 +1,12 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/external_imports.R +\name{\%>\%} +\alias{\%>\%} +\title{Pipe operator} +\usage{ +lhs \%>\% rhs +} +\description{ +See \code{magrittr::\link[magrittr:pipe]{\%>\%}} for details. +} +\keyword{internal} diff --git a/man/threshold_VAR.Rd b/man/threshold_VAR.Rd new file mode 100644 index 0000000..847924f --- /dev/null +++ b/man/threshold_VAR.Rd @@ -0,0 +1,36 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/threshold_var.R +\name{threshold_VAR} +\alias{threshold_VAR} +\title{Estimate multi-regime VAR} +\usage{ +threshold_VAR(data, regime, p = 1, horizon = 10, freq = "month") +} +\arguments{ +\item{data}{data.frame, matrix, ts, xts, zoo: Endogenous regressors} + +\item{regime}{string: name or regime assignment vector in the design matrix (data)} + +\item{p}{int: lags} + +\item{horizon}{int: forecast horizons} + +\item{freq}{string: frequency of data (day, week, month, quarter, year)} +} +\value{ +list of lists, each regime returns its own list with elements \code{data}, \code{model}, \code{forecasts}, \code{residuals} +} +\description{ +Estimate multi-regime VAR +} +\examples{ +\dontrun{ +threshold_VAR( + data = Data, + regime = 'regime', + p = 1, + horizon = 10, + freq = 'month') +} + +} diff --git a/man/threshold_lp_irf.Rd b/man/threshold_lp_irf.Rd new file mode 100644 index 0000000..08d78ff --- /dev/null +++ b/man/threshold_lp_irf.Rd @@ -0,0 +1,39 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/local_projections.R +\name{threshold_lp_irf} +\alias{threshold_lp_irf} +\title{Estimate multi-regime local projection IRFs} +\usage{ +threshold_lp_irf(data, shock, target, thresholdVar = NULL, horizons, lags) +} +\arguments{ +\item{data}{data.frame, matrix, ts, xts, zoo: Endogenous regressors} + +\item{shock}{string: variable to shock} + +\item{target}{string: variable betas to collect} + +\item{thresholdVar}{data.frame of regime binaries} + +\item{horizons}{int: horizons to forecast out to} + +\item{lags}{int: lags to include in regressions} +} +\value{ +data.frame +} +\description{ +Estimate multi-regime local projection IRFs +} +\examples{ +\dontrun{ +threshold_lp_irf( + data = Data, + shock = 'x', + target = 'y', + thresholdVar = Data.regime, + horizons = 20, + lags = 2) +} + +} diff --git a/man/threshold_var_fevd.Rd b/man/threshold_var_fevd.Rd new file mode 100644 index 0000000..b21057c --- /dev/null +++ b/man/threshold_var_fevd.Rd @@ -0,0 +1,27 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/threshold_var_fevd.R +\name{threshold_var_fevd} +\alias{threshold_var_fevd} +\title{Estimate multi-regime forecast error variance decomposition} +\usage{ +threshold_var_fevd(threshold_var, horizon = 10) +} +\arguments{ +\item{threshold_var}{threshold_var output} + +\item{horizon}{int: number of periods} +} +\value{ +list, each regime returns its own long-form data.frame +} +\description{ +Estimate multi-regime forecast error variance decomposition +} +\examples{ +\dontrun{ +threshold_var_fevd( + threshold_var, + bootstraps.num = 10) +} + +} diff --git a/man/threshold_var_irf.Rd b/man/threshold_var_irf.Rd new file mode 100644 index 0000000..8fab27e --- /dev/null +++ b/man/threshold_var_irf.Rd @@ -0,0 +1,37 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/threshold_var_irf.R +\name{threshold_var_irf} +\alias{threshold_var_irf} +\title{Estimate multi-regime impulse response functions} +\usage{ +threshold_var_irf( + threshold_var, + horizon = 10, + bootstraps.num = 100, + CI = c(0.1, 0.9) +) +} +\arguments{ +\item{threshold_var}{threshold_VAR output} + +\item{horizon}{int: number of periods} + +\item{bootstraps.num}{int: number of bootstraps} + +\item{CI}{numeric vector: c(lower ci bound, upper ci bound)} +} +\value{ +list of lists, each regime returns its own list with elements \code{irfs}, \code{ci.lower}, and \code{ci.upper}; all elements are long-form data.frames +} +\description{ +Estimate multi-regime impulse response functions +} +\examples{ +\dontrun{ +threshold_var_irf( + threshold_var, + bootstraps.num = 10, + CI = c(0.05,0.95)) +} + +} diff --git a/man/var_fevd.Rd b/man/var_fevd.Rd new file mode 100644 index 0000000..e88c33e --- /dev/null +++ b/man/var_fevd.Rd @@ -0,0 +1,27 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/var_fevd.R +\name{var_fevd} +\alias{var_fevd} +\title{Estimate single-regime forecast error variance decomposition} +\usage{ +var_fevd(var, horizon = 10) +} +\arguments{ +\item{var}{VAR output} + +\item{horizon}{int: number of periods} +} +\value{ +long-form data.frame +} +\description{ +Estimate single-regime forecast error variance decomposition +} +\examples{ +\dontrun{ +var_fevd( + var, + bootstraps.num = 10) +} + +} diff --git a/man/var_irf.Rd b/man/var_irf.Rd new file mode 100644 index 0000000..3d4e9c6 --- /dev/null +++ b/man/var_irf.Rd @@ -0,0 +1,32 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/var_irf.R +\name{var_irf} +\alias{var_irf} +\title{Estimate single-regime impulse response functions} +\usage{ +var_irf(var, horizon = 10, bootstraps.num = 100, CI = c(0.1, 0.9)) +} +\arguments{ +\item{var}{VAR output} + +\item{horizon}{int: number of periods} + +\item{bootstraps.num}{int: number of bootstraps} + +\item{CI}{numeric vector: c(lower ci bound, upper ci bound)} +} +\value{ +list object with elements \code{irfs}, \code{ci.lower}, and \code{ci.upper}; all elements are long-form data.frames +} +\description{ +Estimate single-regime impulse response functions +} +\examples{ +\dontrun{ +var_irf( + var, + bootstraps.num = 10, + CI = c(0.05,0.95)) +} + +} diff --git a/man/var_irf_plot.Rd b/man/var_irf_plot.Rd new file mode 100644 index 0000000..72e5d0f --- /dev/null +++ b/man/var_irf_plot.Rd @@ -0,0 +1,21 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/var_irf.R +\name{var_irf_plot} +\alias{var_irf_plot} +\title{Plot all IRFs} +\usage{ +var_irf_plot(irfs, shocks, responses) +} +\arguments{ +\item{irfs}{var_irf object} + +\item{shocks}{string vector: shocks to plot} + +\item{responses}{string vector: responses to plot} +} +\value{ +grid of ggplot2 graphs +} +\description{ +Plot all IRFs +} diff --git a/sovereign.Rproj b/sovereign.Rproj new file mode 100644 index 0000000..69fafd4 --- /dev/null +++ b/sovereign.Rproj @@ -0,0 +1,22 @@ +Version: 1.0 + +RestoreWorkspace: No +SaveWorkspace: No +AlwaysSaveHistory: Default + +EnableCodeIndexing: Yes +UseSpacesForTab: Yes +NumSpacesForTab: 2 +Encoding: UTF-8 + +RnwWeave: Sweave +LaTeX: pdfLaTeX + +AutoAppendNewline: Yes +StripTrailingWhitespace: Yes +LineEndingConversion: Posix + +BuildType: Package +PackageUseDevtools: Yes +PackageInstallArgs: --no-multiarch --with-keep.source +PackageRoxygenize: rd,collate,namespace diff --git a/tests/testthat.R b/tests/testthat.R new file mode 100644 index 0000000..9d56026 --- /dev/null +++ b/tests/testthat.R @@ -0,0 +1,4 @@ +library(testthat) +library(sovereign) + +test_check("sovereign") diff --git a/tests/testthat/test-lp.R b/tests/testthat/test-lp.R new file mode 100644 index 0000000..63ac611 --- /dev/null +++ b/tests/testthat/test-lp.R @@ -0,0 +1,33 @@ +test_that("Local projection workflow", { + + # simple time series + AA = c(1:100) + rnorm(100) + BB = c(1:100) + rnorm(100) + CC = AA + BB + rnorm(100) + date = seq.Date(from = as.Date('2000-01-01'), by = 'month', length.out = 100) + Data = data.frame(date = date, AA, BB, CC) + Data = dplyr::mutate(Data, reg = dplyr::if_else(AA > median(AA), 1, 0)) + + # local proejctions + irfs = + lp_irf( + data = Data, + shock = 'AA', + target = 'BB', + horizons = 20, + lags = 2) + + t.irfs = + threshold_lp_irf( + data = dplyr::select(Data, -reg), + thresholdVar = dplyr::select(Data, reg, date), + shock = 'AA', + target = 'BB', + horizons = 20, + lags = 2) + + expect_true(is.data.frame(irfs)) + expect_true(is.list(t.irfs)) + + +}) diff --git a/tests/testthat/test-regimes.R b/tests/testthat/test-regimes.R new file mode 100644 index 0000000..75f9430 --- /dev/null +++ b/tests/testthat/test-regimes.R @@ -0,0 +1,28 @@ +test_that("Regime learning functions", { + + # simple time series + AA = c(1:100) + rnorm(100) + BB = c(1:100) + rnorm(100) + CC = AA + BB + rnorm(100) + date = seq.Date(from = as.Date('2000-01-01'), by = 'month', length.out = 100) + Data = data.frame(date = date, AA, BB, CC) + Data = dplyr::mutate(Data, reg = dplyr::if_else(AA > median(AA), 1, 0)) + + # run ml + regimes.rf = + learn_regimes( + data = Data, + engine = 'rf' + ) + + regime.kmeans = + learn_regimes( + data = Data, + engine = 'kmeans', + regime.n = 2 + ) + + expect_true(is.data.frame(regimes.rf)) + expect_true(is.data.frame(regime.kmeans)) + +}) diff --git a/tests/testthat/test-threvhold_var.R b/tests/testthat/test-threvhold_var.R new file mode 100644 index 0000000..5beb7bf --- /dev/null +++ b/tests/testthat/test-threvhold_var.R @@ -0,0 +1,44 @@ +test_that("threshold VAR workflow", { + + # simple time series + AA = c(1:100) + rnorm(100) + BB = c(1:100) + rnorm(100) + CC = AA + BB + rnorm(100) + date = seq.Date(from = as.Date('2000-01-01'), by = 'month', length.out = 100) + Data = data.frame(date = date, AA, BB, CC) + Data = dplyr::mutate(Data, reg = dplyr::if_else(AA > median(AA), 1, 0)) + + # estimate VAR + tvar = + threshold_VAR( + data = Data, + regime = 'reg', + p = 1, + horizon = 10, + freq = 'month') + + expect_true(is.list(tvar)) + expect_true(is.list(tvar$model)) + expect_true(is.list(tvar$forecasts)) + expect_true(is.list(tvar$residuals)) + + # estimate IRF + irf = + threshold_var_irf( + tvar, + bootstraps.num = 10, + CI = c(0.05,0.95)) + + expect_true(is.list(irf)) + expect_true(is.data.frame(irf[[1]]$ci.lower)) + expect_true(is.data.frame(irf[[1]]$ci.upper)) + + # estimate forecast error variance decomposition + fevd = + threshold_var_fevd( + tvar, + horizon = 10) + + expect_true(is.list(fevd)) + +}) diff --git a/tests/testthat/test-var.R b/tests/testthat/test-var.R new file mode 100644 index 0000000..9533558 --- /dev/null +++ b/tests/testthat/test-var.R @@ -0,0 +1,42 @@ +test_that("VAR workflow", { + + # simple time series + AA = c(1:100) + rnorm(100) + BB = c(1:100) + rnorm(100) + CC = AA + BB + rnorm(100) + date = seq.Date(from = as.Date('2000-01-01'), by = 'month', length.out = 100) + Data = data.frame(date = date, AA, BB, CC) + + # estimate VAR + var = + VAR( + data = Data, + p = 1, + horizon = 10, + freq = 'month') + + expect_true(is.list(var)) + expect_true(is.list(var$model)) + expect_true(is.list(var$forecasts)) + expect_true(is.list(var$residuals)) + + # estimate IRF + irf = + var_irf( + var, + bootstraps.num = 10, + CI = c(0.05,0.95)) + + expect_true(is.list(irf)) + expect_true(is.data.frame(irf$ci.lower)) + expect_true(is.data.frame(irf$ci.upper)) + + # estimate forecast error variance decomposition + fevd = + var_fevd( + var, + horizon = 10) + + expect_true(is.data.frame(fevd)) + +}) diff --git a/vignettes/getting_started.Rmd b/vignettes/getting_started.Rmd new file mode 100644 index 0000000..5a51fa3 --- /dev/null +++ b/vignettes/getting_started.Rmd @@ -0,0 +1,152 @@ +--- +title: "Sovereign, getting started" +output: rmarkdown::html_vignette +vignette: > + %\VignetteIndexEntry{getting_started} + %\VignetteEngine{knitr::rmarkdown} + %\VignetteEncoding{UTF-8} +--- + +```{r, include = FALSE} +knitr::opts_chunk$set( + collapse = TRUE, + comment = "#>" +) + +options(rmarkdown.html_vignette.check_title = FALSE) +``` + + +As the sovereign package is under active development, its **API is not guaranteed to be stable**. As a result, this introductory vignette will remain fairly barebones for the time being, and is meant to simply and succinctly demonstrate the most up-to-date workflow supported by the sovereign package. + +## 0. Environment + +```{r setup} +# load packages +suppressPackageStartupMessages(library(sovereign)) # analysis +suppressPackageStartupMessages(library(quantmod)) # FRED api +suppressPackageStartupMessages(library(dplyr)) # general cleaning +suppressPackageStartupMessages(library(lubridate)) # date functions +``` + +## 1. Create Data + +```{r} +# pull and prepare data from FRED +quantmod::getSymbols.FRED( + c('UNRATE','INDPRO','GS10'), + env = globalenv()) + +Data = cbind(UNRATE, INDPRO, GS10) + +Data = data.frame(Data, date = zoo::index(Data)) %>% + dplyr::filter(lubridate::year(date) >= 1990) %>% + na.omit() + +``` + +## 2. Assign Regimes + +```{r} + +# create a regime explicitly +Data.threshold = Data %>% + mutate(mp = if_else(GS10 > median(GS10), 1, 0)) + +# assign regimes based on unsurpervised kmeans +# (will not be used further) +regimes = + learn_regimes( + data = Data, + regime.n = 3, + engine = 'kmeans') + +``` + +## 3. Single-Regime Analaysis + +```{r} + +#------------------------------------------ +# single-regime var +#------------------------------------------ +# estimate VAR +var = + VAR( + data = Data, + p = 1, + horizon = 10, + freq = 'month') + +# estimate IRF +irf = + var_irf( + var, + bootstraps.num = 10, + CI = c(0.05,0.95)) + +# estimate forecast error variance decomposition +fevd = + var_fevd( + var, + horizon = 10) + +#------------------------------------------ +# single-regime local projection +#------------------------------------------ +# estimate IRF + lp.irfs = + lp_irf( + data = Data, + shock = 'INDPRO', + target = 'GS10', + horizons = 20, + lags = 2) + +``` + +## 4. Multi-Regime Analysis + +```{r} + +#------------------------------------------- +# multi-regime var +#------------------------------------------ +# estimate multi-regime VAR +tvar = + threshold_VAR( + data = Data.threshold, + regime = 'mp', + p = 1, + horizon = 1, + freq = 'month') + +# estimate IRF +tvar.irf = + threshold_var_irf( + tvar, + horizon = 10, + bootstraps.num = 10, + CI = c(0.05,0.95)) + +# estimate forecast error variance decomposition +tvar.fevd = + threshold_var_fevd( + tvar, + horizon = 10) + +#------------------------------------------- +# multi-regime local projection +#------------------------------------------ +# estimate IRF +tlp.irfs = + threshold_lp_irf( + data = dplyr::select(Data.threshold, -mp), + thresholdVar = dplyr::select(Data.threshold, mp, date), + shock = 'INDPRO', + target = 'GS10', + horizons = 20, + lags = 2) + + +```