Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
1,124,944 characters · 357 sections · 115 citation commands
Quantitative methods in finance
\sloppy
{{0.6\textwidth} \epigraph{Perdre un bonheur r\^{e}v\'{e}, renoncer \`{a} tout un avenir, est une souffrance plus aigu\"{e} que celle caus\'{e}e par la ruine d'une f\'{e}licit\'{e} ressentie, quelque compl\`{e}te qu'elle ait \'{e}t\'{e} : l'esp\'{e}rance n'est-elle pas meilleure que le souvenir ?}{La Bourse\\Honor\'{e} de Balzac}} \setcounter{tocdepth}{2}
\quad
\newcounter{question}
This course, Quantitative Methods in Finance, is designed for Master 2 Research students and provides a comprehensive introduction to advanced econometric and statistical methods with direct applications in finance. The objective is to equip students with practical quantitative skills relevant for careers across the financial sector, including banking and financial services, asset management, fintech, start-ups and artificial intelligence, (re)insurance, and public financial institutions and central banking.
We combine theoretical foundations with hands-on Python implementation, introducing programming from scratch and covering modern techniques such as large language model (LLM) prompt engineering for financial data analysis. Students will learn to collect and process data from traditional (central banks, regulators, international organisations, stock exchanges) and alternative sources, and to implement and interpret advanced econometric models.
Topics include: data handling and quality assessment; probability theory and statistical inference; regression analysis; time series modelling (stationarity, seasonality, SARIMAX, GARCH); structural VARs and impulse–response analysis; Difference-in-Differences and event studies; risk modelling via Extreme Value Theory; and dependence modelling with copulas. A dedicated session situates econometrics within the broader machine learning landscape, highlighting when to prefer parametric approaches over more flexible ML methods.
By the end of the course, students will be able to replicate state-of-the-art empirical papers and carry out their own econometric analyses. Assessment is based on in-class exercises and a final project.
Two open source tools are suggested to be used during the class: python and R. While both tools have pros and cons (\href{http://www.tiobe.com/tiobe-index/}{and their rankings depends on which criteria you value, e.g. based on search engines}), they both allow for rapid prototyping of investment strategies and/or academic models which is our main objective. In the industry, once back-tested an investment strategy would be implemented typically in C++.
Students are also allowed to use Julia, Matlab, Stata or SAS for their work as long as they can demonstrate on request they have a valid license if the software requires one. If a pirate software is found on a computer, sanctions will apply.
We strongly recommend students to invest in a laptop. As of 2026, based on the website \href{https://www.lesnumeriques.com/ordinateur-portable/guide-achat-quels-sont-les-meilleurs-ordinateurs-portables-g251.html}{lesnumeriques} it is feasible to have a computer well performing for a budget around 800 \euro, when choosing a laptop keep in mind that rapid access to the hard drive is recommended (which favour SSD), if one wish to work with big data the RAM will be important (16 Go is a good starting point), then for computation the CPU enters into consideration (1.5 GHz is a good starting point).
As for the LLM question, there exist benchmarks. As of Sep. 2025, I am using the "Plus" version at 23 euros per month. In the year 2024-2025, a student had access to the Pro 229 euros per month version. This challenges us, teachers, to be able to evaluate the student's contribution controlling for the LLM version used. This is pointing toward more in-class pen and paper tests, as in a recruitment interview for a quantitative role.
The aim of these sessions is not (yet) to produce innovations, but rather to learn to replicate academic research papers using available data (students should register at the library where they would have access to data sources). Once the quantitative methods are mastered, the students should be able to contribute to the economic research field. The students should be able to solve basic quantitative exercises with pen and papers.
Eric Vansteenberghe is a research economist with more than ten years of experience teaching quantitative methods to Master 2 research students. His research agenda focuses on financial intermediation, monetary policy, uncertainty, and risk, with applications spanning banking, credit markets, financial technology, and insurance. A central contribution of his work is the empirical identification of mechanisms operating under uncertainty and asymmetry, using granular data and transparent econometric designs. This includes analyses of how monetary policy interacts with uncertainty to shape credit supply vansteenberghe2025monetary and vansteenberghe2024uncertain, how public guarantees affect bank risk-taking and regulatory capital incentives nicolas2025public, and how financial technologies modify collateral use and lending outcomes for small businesses beaumont2025collateral, in close connection with the broader literature on FinTech and banking competition doerr2024literature. His research also addresses insurance supervision and climate-related risks, proposing new detection methods to identify vulnerabilities arising from climate change within insurance systems vansteenberghe2025insurance. The Quantitative Methods in Finance (QMF) lecture notes are the direct outcome of this research–teaching nexus. They consolidate more than a decade of material developed for research-oriented Master students and are deliberately anchored in the author’s published work and methodological standards. In doing so, the notes aim not only to serve as a pedagogical reference but also to provide large language models and other automated tools with explicit academic grounding, references, and conventions consistent with contemporary research practice.
The methods and conclusions reached in this documents and codes are:
According to the Financial Times of 3rd February 2020, Citigroup has suspended one of its most senior bond traders in London (Paras Shah) after he had stolen food from the canteen. Two of his former colleagues told the FT he was a well-liked and successful trader.
In 2016, Japan's Mizuho Bank fired a London banker after he was caught stealing a part from a colleague's bike worth about 5 GBP.
I first found it hard to understand how wealthy bankers would steal a sandwich or a bicycle bell. But then I realized how, in a rush, we might misjudge a situation and be tempted to take short-cuts. It would have taken both employees half an hour to go to a regular shop to buy a sandwich or a bicycle bell. They could both afford spending half an hour on this task, or even 5 minutes to order it online, and spending the money. They did not bother and they lost their jobs and their names are out there, making it hard for them to find another position at a bank.
This is a reminder to be humble and even when it is 20:30 and we are tired not to take shortcuts that can cost our reputation.
The students will be evaluated via several methods. The final grade is out of 20.
Five points will be given to in-class pen and paper exercises related to the quantitative work of this lecture.
Bonus points can be given based on in-class exercises. Exercises have to be carried out by the students and can be requested at the end of the lecture to be evaluated. Alternatively, the students can be selected to present a technical subject and programming application during the lecture. If she has not been selected during the semester at the end of a class or for a presentation, the students can have to work during the final session on coding in R or python. Work will be individual. Internet connexion might not be permitted during the whole session (to be determined on the day depending on the subject). In case the student has no computer(or its battery is empty) or its Spyder or RStudio environments are not working she will have to submit her work in paper form, otherwise submission is expected at the end in e-mail format to the supervisor.
Fifteen points are allocated to the final report. For the final report, students have to prepare a summary of the theory (for example re-demonstrate some key equations to clarify how the model work and on which theory it is based), and to reproduce the models and results of research papers to address a question (for this semester, the same question is given to all students and pairs). Each student will have to:
To be announced.
We want to test the following skills:
If you want a reference on how to best write your code in python, you can follow \href{https://www.python.org/dev/peps/pep-0008/}{PEP 8}.
As a minimum, in your python or R codes:
Nota bene:
This section is complemented by original teaching material developed for this course by Eric Vansteenberghe and made publicly available for replication and self-study. In particular, the GitHub repository \href{https://github.com/skimeur/QMF}{github.com/skimeur/QMF} contains the full set of course resources described in the README: structured lecture notes, Python and R scripts, reproducible empirical exercises, environment setup instructions, and worked examples designed to bridge economic theory, econometrics, and applied data analysis. The repository is intended as a living companion to the lectures, allowing students to reproduce results, extend exercises, and progressively build autonomous research workflows.
The theoretical foundations and pedagogical structure of the course are documented in the lecture notes vansteenberghe2026quantitative, which are available both as an \href{https://arxiv.org/abs/2601.12896}{arXiv preprint} and on \href{https://papers.ssrn.com/sol3/papers.cfm?abstract_id=5178205}{SSRN}. These notes provide a unified treatment of quantitative methods in finance, with an emphasis on empirical implementation in Python and R, and serve as the primary written reference for the material covered in this section.
If you have never used any tool like R, Matlab, Stata, you might want to start with those step by step exercises we suggest to understand how to use matrices, DataFrames, functions, loops:
\quad
The aim of our Quantitative Methods in Finance is to provide students with tools to deploy the economic and mathematical theories to empirical data sets. There are some resources if students want to learn the languages R and python in depth:
The course material is hosted on GitHub. The expected workflow is the standard open-source contribution cycle: fork the repository, implement your changes on your fork, and propose them back to the upstream repository so that I can review, request edits if needed, and merge your contribution. Once merged, your work becomes part of the official codebase, visible to future cohorts (with attribution through Git history and the pull request).
\paragraph{Step 1 --- Fork the repository (create your own copy).}
\paragraph{Step 2 --- Create a clean working branch.} Do not work directly on main. Create a dedicated branch named after your task:
\paragraph{Step 3 --- Implement your changes (reproducible and documented).}
\paragraph{Step 4 --- Commit your work with clear messages.} Commit early and often. Use informative commit messages:
\paragraph{Step 5 --- Push your branch to GitHub.}
\paragraph{Step 6 --- Open a Pull Request (PR) to propose your changes.} Once you finished all your work and are ready to submit:
\paragraph{Step 7 --- Review cycle (what happens after you open the PR).}
\paragraph{Attribution and long-run value.} After merge, your contribution is permanently visible:
\paragraph{Good practice checklist (to maximize acceptance probability).}
In this course, we will use a public repository maintained by the instructor:
Students can read and replicate everything, but do not have direct write access. Contributions therefore follow the standard open-source workflow: fork $\rightarrow$ modify $\rightarrow$ pull request.
\paragraph{Repository and Local Copy.} A repository is a project space that contains all files together with their complete revision history. In Git, repositories are distributed: each local clone contains the full history, while GitHub hosts a public reference version (the upstream repository). Before cloning a repository, you are free to choose where on your computer the project will be stored (folder and path).
\paragraph{Checking the Current Path.} To display the current working directory in the terminal, use:
This command (print working directory) shows the full path of the folder in which you are currently located.
\paragraph{Choosing the Local Folder.} First, open a terminal and navigate to the directory in which you want the project to be stored:
For example:
\paragraph{Contribute Improvements (Fork-and-Pull Workflow).} If you want to propose an improvement that may be integrated into the reference repository, you must use a fork and submit a pull request. This keeps all contributions public and ensures the instructor retains editorial control.
\paragraph{Step 1: Fork the Repository on GitHub.} On GitHub:
\paragraph{Step 2: Clone Your Fork Locally.} Clone the repository from your own GitHub account (not the instructor's):
\paragraph{Step 3: Create a Dedicated Branch.} Create a branch describing your contribution (recommended):
Using branches keeps your work isolated and makes review easier.
\paragraph{Step 4: Edit Files and Check the Project Status.} Modify files locally, then inspect what changed:
To inspect the exact changes:
\paragraph{Step 5: Stage and Record Changes (Commit).} Stage your modifications:
or stage everything you intentionally changed:
Then commit with an informative message:
\paragraph{Step 6: Push Your Branch to Your Fork.} Upload your branch to your GitHub fork:
\paragraph{Step 7: Open a Pull Request (PR).} On GitHub:
\paragraph{Step 8: Review and Decision (Instructor Side).} The instructor will review your pull request and may:
Only pull requests explicitly merged by the instructor modify the reference repository used by the class.
\paragraph{Synchronizing Your Fork with the Instructor Repository (Optional).} Over time, the instructor repository will evolve. To keep your fork up to date, add the instructor repository as an upstream remote once:
Then periodically fetch and integrate updates from the instructor repository:
This operation updates your fork without modifying the instructor repository.
Remark. Advanced users may prefer git rebase upstream/main instead of merge to maintain a linear history, but this is not required for the course.
\paragraph{Good Practices for Student Contributions.} To ensure reproducibility and clarity, contributors should:
Mastering these steps is sufficient to (i) replicate the repository locally and (ii) propose public improvements under full editorial control of the maintainer.
Version control systems allow researchers to track changes in code, data, and manuscripts over time, even when working alone. Git records changes locally, while GitHub provides a remote repository that enables backup, synchronization across machines, and future collaboration.
In this subsection, GitHub is used as a personal research archive: a repository controlled by a single user, either private or public.
\paragraph{Creating a Repository on GitHub.} Create a new repository from your GitHub account:
GitHub will provide the repository URL for linking your local project.
\paragraph{Repository and Local Copy.} A repository is a structured folder containing project files and their complete revision history. You are free to choose the location of the local project folder on your machine.
\paragraph{Checking the Current Path.} To display your current working directory:
\paragraph{Choosing the Local Folder.} Navigate to the directory where you store research projects:
For example:
\paragraph{Option A: Cloning an Existing Repository.} If the repository was created on GitHub first:
This creates a local folder already linked to GitHub.
\paragraph{Option B: Initializing Git in an Existing Folder.} If the project already exists locally:
\paragraph{Continuing Work on an Already Cloned Project.}
When reopening an existing project, always ensure that your local copy is synchronized with GitHub before starting new modifications:
Step 1: Switch to the main branch.
Step 2: Pull the latest changes.
This step is essential if:
Failing to pull before editing may create avoidable merge conflicts.
\paragraph{Checking the Project Status.} Before recording changes:
This shows modified, staged, and untracked files.
\paragraph{Staging and Recording Changes.} Stage modifications:
Record a snapshot:
Each commit represents a documented and reproducible step.
\paragraph{Uploading Changes to GitHub.} Push committed changes:
If pushing a branch for the first time:
GitHub then serves as both backup and synchronization device.
\paragraph{Inspecting Past Modifications.} To inspect commit history:
To inspect differences:
\paragraph{Typical Individual Workflow.}
\paragraph{Good Practices for Individual Research.}
Used in this disciplined way, GitHub provides traceability, reproducibility, and a seamless transition to collaborative workflows.
When collaborators have direct write access to the same repository (e.g.\ research team, co-authors), the workflow differs from the fork-and-pull model. The objective is to avoid conflicts, preserve a clean history, and ensure reproducibility.
\paragraph{Recommended Workflow (Shared Repository).}
\paragraph{Rules for Safe Collaboration.}
This workflow minimizes merge conflicts, preserves a transparent research history, and ensures that collaborative academic projects remain reproducible and auditable.
When only two researchers collaborate closely and communicate continuously, working directly on main can be acceptable. This setup is common in small academic projects where both contributors understand the full codebase and coordinate changes in advance.
\paragraph{Minimal Safe Workflow (Direct-to-main).}
\paragraph{Coordination Rules.}
\paragraph{When to Revert to Branching.} If the project grows (teaching use, external replication, additional co-authors, major restructuring), switch to a branch-based workflow to protect stability and traceability.
For small, well-synchronized research teams, disciplined direct commits to main can remain efficient while preserving reproducibility.
This subsection provides a conceptual roadmap of the tools introduced in the following pages. Its purpose is to clarify what each component is, why it is needed, and how the pieces fit together in a modern academic workflow. The perspective adopted is that of a student or researcher with no prior experience in Python or software development, but who needs to (i) write code, (ii) analyze data, and (iii) collaborate reproducibly with others.
\paragraph{Python as the Core Language.} Python is a general-purpose programming language that has become a standard in data analysis, econometrics, machine learning, and scientific computing. It is used to write scripts that load data, perform computations, estimate models, and generate figures and tables. For beginners, Python is attractive because of its readable syntax and extensive documentation. An official entry point is the Python Beginners Guide:
\paragraph{Why Package and Environment Management Is Necessary.} Python itself is only the language. Most practical work relies on external packages (for example, numpy, pandas, matplotlib, statsmodels). These packages evolve over time and may depend on specific versions of other libraries. Without a management system, installing or updating packages can easily break existing code. This is particularly problematic in academic settings, where reproducibility across students, coauthors, and over time is essential.
\paragraph{Conda: Managing Python and Its Dependencies.} Conda is both a package manager and an environment manager. It addresses two key problems:
Each Conda environment can have its own Python version and its own set of packages, ensuring that different projects do not interfere with each other. The official Conda documentation provides a structured introduction:
\paragraph{Anaconda as a Turnkey Distribution.} Anaconda is a complete Python distribution built around Conda. It bundles:
For beginners, Anaconda reduces installation complexity and avoids many configuration pitfalls. The official getting-started guide is available at:
\paragraph{Integrated Development Environments (IDEs).} Writing and maintaining code is done most efficiently within an Integrated Development Environment. An IDE provides syntax highlighting, error detection, code execution, and debugging tools.
Two IDEs are particularly relevant in this course:
Spyder installation instructions are available at:
VS Code documentation is available at:
\paragraph{Interactive Computing with Jupyter.} In addition to scripts and IDEs, Python supports interactive workflows through Jupyter Notebooks. Notebooks allow code, results, figures, and explanatory text to coexist in a single document. They are particularly useful for:
Modern IDEs (Spyder and VS Code) integrate Jupyter-style execution directly into standard .py files or notebook interfaces.
\paragraph{Version Control and Collaboration with Git and GitHub.} When working with others, code and data must be shared, tracked, and synchronized. This is handled through:
GitHub is not only a storage service: it is a central tool for collaborative academic research, allowing multiple contributors to work on the same codebase while preserving a complete history of modifications.
\paragraph{Putting Everything Together: A Typical Workflow.} In practice, a standard academic workflow combines all these components:
The following subsections provide step-by-step instructions for installing and using each of these tools. The goal is not only to make the code run, but to establish good research practices in terms of reproducibility, clarity, and collaboration.
For this course, we recommend relying on Conda for Python distribution, package management, and—most importantly—environment management. Conda provides a controlled and reproducible way to install Python and its scientific ecosystem without interfering with the system Python installation.
\paragraph{What is Conda?} Conda is an open-source package manager and environment manager. It serves two distinct but related purposes:
This separation is crucial in practice\footnote{Although to be transparent, I never had to worry about it as my research projects are short enough not to have environment/replicability over time issues.}. Different projects often require incompatible package versions; Conda environments allow these projects to coexist cleanly on the same machine.
\paragraph{Anaconda vs. Miniconda} Anaconda is a full distribution built on top of Conda. It ships with Python, Conda, and a large collection of pre-installed scientific packages (NumPy, pandas, SciPy, matplotlib, etc.), as well as graphical tools such as Anaconda Navigator. While lighter alternatives exist (e.g.\ Miniconda), Anaconda is the preferred solution for this course because it minimizes setup frictions and provides an integrated visual interface.
Install Anaconda for Python 3 from the official website: \[ \texttt{\href{https://www.anaconda.com/download}{https://www.anaconda.com/download}} \] This installation includes:
No separate Python installation is required.
A Conda environment is a self-contained directory that holds:
Environments are independent: installing or updating a package in one environment has no effect on the others. This is essential for reproducibility and for avoiding dependency conflicts.
Typical workflow:
Packages are installed from channels (repositories), with Conda automatically resolving dependencies. Compared to standard Python tools, Conda can also manage non-Python dependencies (compiled libraries), which is particularly important for scientific computing.
Conceptually:
For a guided introduction to the IDE, see Help $\rightarrow$ Spyder tutorial.
Two complementary tools are commonly used when working with Python: Integrated Development Environments (IDEs) and Jupyter Notebooks. They serve different purposes and are often used jointly within the same project.
\paragraph{IDEs (e.g.\ VS Code).} An IDE such as Visual Studio Code is designed for structured software development. According to widely used tutorials and official documentation, its main features include:
IDEs are particularly well suited for writing reusable scripts, managing multi-file projects, and ensuring reproducibility in research and production settings.
\paragraph{Jupyter Notebooks (e.g.\ Google Colab).} Jupyter Notebooks provide an interactive, cell-based execution model that combines code, outputs, figures, and explanatory text. As emphasized in most introductory tutorials, they are especially useful for:
Cloud-based platforms such as Google Colab allow users to run notebooks without local installation, offering immediate access to a preconfigured environment, at the cost of more limited control over package versions and long-term reproducibility.
\paragraph{Choosing between IDEs and Notebooks.} In practice, notebooks are typically favored for exploration and communication, while IDEs are preferred for developing, structuring, and maintaining code. Modern workflows often combine both approaches: ideas are prototyped in notebooks and later translated into scripts or modules managed within an IDE.
To work efficiently with GitHub repositories, we recommend using Visual Studio Code (VS Code), a lightweight and widely adopted code editor that integrates seamlessly with Git, Python, and \LaTeX.
\paragraph{Installing Visual Studio Code.} Visual Studio Code is freely available for all major operating systems (Windows, macOS, Linux). It can be downloaded from the official website:
After installation, VS Code can be launched as a standard application.
\paragraph{Installing required VS Code extensions (first use).} Before running Python code interactively, two extensions must be installed in VS Code:
These extensions can be installed via the Extensions panel in VS Code. The Python extension enables script execution and interpreter selection, while the Jupyter extension is required for interactive execution, cell-based workflows, and variable inspection (even when working with standard .py files).
\paragraph{Opening a Terminal in VS Code.} VS Code provides an integrated terminal, allowing users to execute Git commands directly within the editor. The terminal can be opened via:
or using the shortcut Ctrl+` (Windows/Linux) or Cmd+` (macOS).
Visual Studio Code provides an integrated terminal that allows you to run standard command-line tools directly inside the editor. Git commands can be executed from this terminal in the same way as from a system terminal.
\paragraph{Checking Whether Git Is Installed.} Before cloning a repository, verify that Git is available on your system:
If a version number is displayed, Git is already installed and ready to use.
\paragraph{Installing Git (If Needed).} If the command is not recognized, Git must be installed:
\paragraph{Downloading a GitHub Repository Locally.} To obtain a local copy of an existing GitHub project, users must first choose a directory on their computer where the project will be stored. In the integrated terminal, navigate to this directory and run:
This command downloads the full repository, including all files and version history, into a new folder named pioneer-detection-method.
When a repository is downloaded using git clone, it is saved in the directory where the terminal was located at the time the command was executed. To display this location, use:
The project folder pioneer-detection-method is located inside the directory printed by this command.
\paragraph{Opening the Project in VS Code.} Once the repository has been cloned, the project can be opened in VS Code by selecting:
and choosing the newly created directory. All source files, scripts, and documentation are then accessible from the file explorer panel.
\paragraph{First steps in VS Code after opening the repository.} Once the folder pioneer-detection-method is opened in VS Code, the following steps are required before running the code.
If the environment is correctly configured, the script will execute and generate the corresponding output files (PDF figures).
\paragraph{Running code step-by-step and inspecting variables (VS Code).} To reproduce an interactive workflow similar to Spyder (e.g.\ executing lines with F9 and inspecting variables), code must be executed in the Python Interactive Window, which relies on the Jupyter extension.
\noindentImportant: If Shift+Enter does not open the Python Interactive Window or variables are not visible, verify that both the Python and Jupyter extensions are installed and enabled in VS Code.
Some GitHub repositories store large data files using Git Large File Storage (Git LFS) rather than standard Git tracking. In this case, the files downloaded at cloning time may only contain a small text pointer instead of the actual data.
\paragraph{Typical Symptom.} If a data file opens as a short text such as:
then the real file content has not been downloaded yet.
\paragraph{Required Setup.} Before working with such repositories, install and activate Git LFS:
(installation is required only once per computer).
\paragraph{Downloading the Actual Data.} After cloning the repository, run from the project root:
This command replaces pointer files with the actual large data files.
\paragraph{Important Note.} Without Git LFS properly installed, code may fail silently or produce cryptic errors, since programs (R, Python) will read pointer files instead of real datasets. Students should therefore always check whether a project relies on Git LFS when working with large files.
You already have R and RStudio installed on your computer. This is sufficient: Visual Studio Code will simply reuse your existing R installation.
\paragraph{Install Visual Studio Code Extensions.} Open Visual Studio Code and install the following extensions:
Visual Studio Code itself is developed by :contentReference[oaicite:0]{index=0}.
\paragraph{Link Visual Studio Code to Your R Installation.} Since R is already installed, VS Code should detect it automatically. If not, specify the R executable path in the VS Code settings:
(adapt the path if R is installed elsewhere).
\paragraph{Running R Code.} Open an .R file. You can:
The code is executed in an R terminal embedded directly in Visual Studio Code.
\paragraph{Installing R Packages.} Packages are installed exactly as in RStudio, from the R console:
To install several packages at once:
Installed packages are shared across RStudio and Visual Studio Code, since both rely on the same R installation.
\paragraph{Project-Based Work (Recommended).} If the GitHub repository contains an R project (.Rproj) or uses a package manager (e.g.\ renv), open the project folder in Visual Studio Code and run:
This ensures package versions are consistent across collaborators.
\paragraph{Key Point.} Visual Studio Code does not replace R or RStudio. It is an alternative interface that allows you to work on R code within a broader development environment, especially convenient when collaborating on GitHub projects mixing R, Python, and \LaTeX.
For the purpose of the following lecture section, we recommend to install the latest version of R and a suite that allows to leverage on the power of many packages and offers a visual environment: RStudio Desktop (\href{https://rstudio.com/}{click here to be redirected}). NB: although it is technically possible to run RStudio from Anaconda, we do not recommend this as we have experienced issues in the past and recommend to follow the steps below to install R and RStudio as a standalone.
Install the latest R version from \url{https://cran.r-project.org/}
Install RStudio Desktop
Open RStudio
For a full review of the types in python, you can go to the documentation: \href{https://docs.python.org/3/library/stdtypes.html}{here for the built-in types and the operations you can apply to them.}\footnote{No need to learn the operations by heart, if you want to do an operation, go and type into Google: \href{https://stackoverflow.com/questions/4344017/how-can-i-get-the-concatenation-of-two-lists-in-python-without-modifying-either}{"python concatenate two lists"}}
We intend to present what you are most likely to manipulate as an economist\footnote{vansteenberghe_types_loops_functions.py}.
Imagine you want to define a variable with name "growth". This variable could be:
If at any point in time you want to clean the "Variable explorer" (and the variables in memory), type in the IPython console %reset and confirm typing y.
If you want to compute growth from a time series, you won't just manipulate multiple variables, you can use mutable build-in python types:
We might want to apply an action if a condition is met or loop a variable through different values or while a condition is met.
We can define our own functions and call that function any time we need to do the same operations as defined in the function and return the results.
You can also chose to work with libraries: numpy or pandas which are powerful and each has its advantages. Just as an illustration we compare the speed to perform matrix multiplication with numpy or pandas objects.
As an economist, you might not need to know how python manage immutable and mutable types variables and memory, \href{https://medium.com/@tyastropheus/tricky-python-i-memory-management-for-mutable-immutable-objects-21507d1e5b95}{here are some elements} and \href{http://foobarnbaz.com/2012/07/08/understanding-python-variables/}{a comparison with C}.
The main takeaway is aliasing where two variable refer to the same object, this link can be broken with .copy(). There are some subtleties, however, which is important to understand for some application in pandas and the management of \href{https://stackoverflow.com/questions/15667189/what-is-the-id-function-used-for/15667328}{object identity} as revealed by some limitation in apply .copy(deep=True) to a DataFrame as detailed \href{https://stackoverflow.com/questions/46327494/python-pandas-dataframe-copydeep-false-vs-copydeep-true-vs}{here}.
But do not worry too much about this, I never encountered this problem and I believe that for your master work you are unlikely to encounter such issues.
We suggest a code\footnote{numpy_exercise_part1.py} to get familiar with some numpy basics features in python.
We suggest in this lecture to work with the spyder environment.
First we want to import the package numpy.
With the package numpy you already have a lot of functions implemented and ready to use.
For example, the function exponential: $$\mbox{exp} : x \to \exp(x)$$ read "funtion exponential that to a variable x associates the value "exponential of x".
If we want to know the value of $exp(1.5)$:
Now we want to use x as a vector containing different values and then for each value of x compute the associated y=f(x), with f being a function (in our case we work with exponential function).
We create a vector called x, ranging from 0 to 9 with steps of 1, either we define manually the array or define a range.
This creates a vector:
Then we do $y = f(x)$, the function $f$ in our case being the exponential.
This returns a vector as output:
Now we want to use a package to be able to plot $y=f(x)$.
We might want a smoother plot, in this case, we will use more computing power and have smaller steps between $0$ and $10$, as an illustration we want to plot by steps of $0.1$.
Now you should create your own function, that you could call $\mbox{my\_f}$: $$\mbox{my\_f} : x \to x^2+x-2$$
Once defined, you want to apply your function to the original $x$ vector to get an output vector $y_f$ and plot it.
Now we want to solve for $f(x)=0$. We can try if a=b using the operator "==":
we find that this is not the case: FALSE
we find that this is the case: TRUE
Indeed, $x^2+x-2=0$ when $x=1$.
Now we can loop through all the element in the vector x ($0,1,\ldots,9$) and print which element solve: f(element)=0.
The good news is that there is already a function to do the above: fsolve.
You need to first import the function from the right package. Now to find the solution to $f(x)=0$, it is in one line:
Not that ",0)" here tells fsolve where to start searching for the root in the abscisse line.
You can get help on this function using the console:
x0 = 0 is the starting estimate for the roots of 'func(x) = 0'
In fsolve documentation, you'll see that fsolve uses MINPACK which itself uses Powell's dog leg method. I am afraid the documentation is quite long, if interested, \href{https://cds.cern.ch/record/126569/files/CM-P00068642.pdfhttps://cds.cern.ch/record/126569/files/CM-P00068642.pdf}{here}.
Define a function: $f(x)=\cos(x^2)$, plot it and find the solution of $f(x)=0$ close to $x=10$.
The idea of a fixed point is simple: a number $x$ is called a fixed point of a function $f$ if it satisfies \[ f(x) = x. \] To approximate such a point, we can use an iterative procedure. Starting from an initial guess $x_0$, we repeatedly apply $f$ to generate a sequence:
This process is called fixed-point iteration. If the sequence $(x_n)$ converges, then its limit $x$ must be a fixed point of $f$.
A famous historical application is the computation of square roots. Suppose we want $\sqrt{a}$ for some $a > 0$. Consider the function
If we start from an initial guess $x_0 > 0$ and apply the iteration \[ x_{n+1} = f(x_n), \] the sequence rapidly converges to $x = \sqrt{a}$, no matter which positive $x_0$ we choose. This procedure is known as the Babylonian method.
The Babylonian method is actually a special case of Newton’s method for solving nonlinear equations. Suppose we want to find a root of a differentiable function $f$, i.e.\ solve $f(x) = 0$. Newton’s method defines the iteration
If the sequence converges to $x$, then $f(x) = 0$.
To see the link, take $f(x) = x^2 - a$. The Newton update (ref) then gives exactly the Babylonian iteration above. This shows how fixed-point iteration and Newton’s method are closely related tools for finding solutions.
The Banach fixed-point theorem (also called the contraction mapping theorem) is a fundamental result in analysis. It states:
\paragraph{Example: iteration with the cosine function.} Consider the function \[ f(x) = \cos(x). \] It is a contraction on the real line around the interval $[0,1]$ since $|f'(x)| = |\!-\sin(x)| \leq 1$, and in fact is strictly less than $1$ near the fixed point.
Thus, if we apply the iteration \[ x_{n+1} = \cos(x_n), \] starting from any initial guess $x_0$, the sequence $(x_n)$ will converge to the unique fixed point of $\cos(x)$, which is approximately $x^* \approx 0.739$.
\paragraph{Error bound.} One of the powerful consequences of Banach’s theorem is that it does not only guarantee convergence, but also gives a bound on the error. If $q$ is a contraction constant (a bound on $|f'(x)|$), then after $n$ steps \[ |x_n - x^*| \;\leq\; \frac{q^n}{1-q}\, |x_1 - x_0|. \]
For instance, if we start with $x_0=1$, the effective contraction factor is about $q \approx 0.85$. This means the error shrinks at least as fast as $0.85^n$, giving rapid convergence to the fixed point.
\stepcounter{question}
If you want to find the minimum of a function $f(x)$, always move in the direction that decreases the function value the fastest. This direction is given by the negative gradient of the function.
\paragraph{Step-by-step intuition.}
\paragraph{Update rule.} Mathematically, this gives the iterative formula: \[ x_{n+1} \;=\; x_n - \gamma \, \nabla f(x_n). \] Here,
\paragraph{Why it works.} The gradient tells us the direction of steepest increase of the function. By subtracting it, we move in the direction of steepest decrease. Repeating this procedure makes the sequence $\{x_n\}$ move closer and closer to a local minimum of $f(x)$.
\paragraph{Key insight for practice.}
As an illustration, we use the gradient descent to find the minimum of the function $f(x) = x^4 - 3x^3 +2$.
We suggest a code\footnote{numpy_exercise_part2.py} with further exercises to get familiar with python's numpy library, we will approximate functions, their derivatives and integrals.
One last visual application: we might want to convince ourselves that: $$ \exp(x) = \lim_{n \to \infty} \left( \left(1 + \frac{x}{n} \right)^n \right)$$
For this, we can simply plot the exponential function against $ \left(1 + \frac{x}{n} \right)^n $ for different values of $n$, for example $[10,100,1000]$.
Then finally, we can convince ourselves visually that: $$ \exp( i \pi) = -1$$
We draw an ellipse:
The equation of such a standard ellipse centred at the origin, with height 2b and width 2a is:
It can be shown that the area inside the ellipse is: $\pi a b$
As an exercise we will compute an approximate of this area by integral. We reorganize equation (ref):
We can compute the area of the semi-ellipse (above part) by:
We want to approximate integral (ref) on the interval $[-a,a]$. We can use the lower and upper Darboux sums which is a simple approach in the spirit of the Riemann integral:
Those two sums, as $y(x)$ is continuous, converge to the integral we want to compute.
\stepcounter{question}
A Generalized Mittag-Leffler function can be defined as:
with $\alpha,\beta \in \mathbb{C}, \mathcal{R}(\alpha)>0, \mathcal{R}(\beta)>0, z \in \mathbb{C}$ and the function Gamma: $\Gamma(r) = \int_0^{+\infty}u^{r-1}e^{-u}du$.
Equation (ref) can be set to approximate the exponential function: $$ E_{1,1}(z)=\sum_{k=0}^\infty \frac{z^k}{k!} $$
We have in practice to replace the sum to a limit $nlim$, instead of $\infty$, that we set to $10$ (red line) and $15$ (dashed red line) for illustration purpose:
For example, Mittag-Leffler type functions combined with Laplace transform are useful to find explicit expressions of insurance companies ruin probability as in Constantinescu2018aa
Remember that the derivative $f'(x)$ of a function $f()$ at $x$ is in fact a rate of change defined as: $$f'(x) = \lim_{h\to 0} \frac{f(x+h)-f(x)}{h}$$
We can convince ourselves that $\ln'(x) = \frac{1}{x}$
We can set our own function to approximate the derivative and check for decreasing values of $h$ if indeed our results are fitting with the function $\frac{1}{x}$, which as expected is the case.
We might want to test with python that: $$\int_{1}^{2} \ln(x) dx = \ln(2) - \ln(1) = \ln(2)$$
And we can approximate: $$\int_{1}^{2} \frac{1}{x} dx = \sum_{i=0}^{N} \gamma (1 + h_i)^{-1}$$
with $\forall i > 0, h_i - h_{i-1} = \gamma $ and $h_N = 1$.
Finally, we can approximate the sinus function with an integral of a Mittag-Leffler function:
$$sin(x) = \int_0^x E_{2,1}(-s^2)ds$$
To decompose, in black we draw the $E_{2,1}(x)$ function. To compute $sin(x_i)$ we create intermediary points (3 here) from 0 to $x_i$ and we add (or subract when negative) the area of the rectangles. This then gives us $sin(x_i)$ represented by a green dot:
This has to be done for every $x_i$ between $0$ and $14$ with a chosen step number (the more steps the more precise but the longer it will take to compute).
\stepcounter{question}
In 1881, Joseph Jagger, an English businessman in the textile sector, visited the Casino of Monte Carlo. Having observed that spinning wheels were never perfectly balanced, Jagger wondered about the potential biases of the casino's roulette wheels. Using empirical averages from his observations and relying on the law of large numbers, he detected significant biases in some roulettes after a month. This insight, coupled with the central limit theorem, allowed him to calculate the estimation error and recognize that certain numbers had a higher probability of winning. Leveraging these findings, he won 2 million francs, equivalent to 10 million euros today.
The Monte Carlo method originates from this empirical approach. It is a tool to estimate averages (or expectations) by simulations, especially when direct calculations are complex or unfeasible. Its application became more prominent with the advent of computers, notably for simulating neutron collisions post-war. Nicholas Metropolis christened this technique in reference to the casino.
At its core, the Monte Carlo method is a statistical simulation technique that leverages the power of random sampling to approximate numerical results. Here are the main three steps:
In the world of econometrics, the Monte Carlo method finds application in various ways. It can be used to assess the properties of an estimator under different conditions, test the power of a statistical test, or even to approximate complex integrals in likelihood functions.
To illustrate\footnote{numpy_exercise_part3.py} how empirical sampling approaches theoretical distributions, we consider the standard normal distribution $N(0,1)$ with probability density function \[ f(x)= \frac{1}{\sqrt{2\pi}} e^{-\tfrac{x^2}{2}}. \] Using random draws of increasing size $n \in \{10^2,10^3,10^4,10^5\}$, we compare the empirical histograms with the true PDF. As $n$ grows, the empirical histogram converges toward the smooth bell curve of the standard normal, highlighting the law of large numbers in action.
In addition to the visual comparison, we compute the probability of observing a draw below $-2$. The theoretical value is given by the standard normal cumulative distribution function, \[ \Pr[X < -2] = \Phi(-2) \approx 0.0228. \] For each sample size, we calculate the empirical frequency of values below $-2$ and report the absolute error relative to the theoretical probability. As expected, the empirical estimates fluctuate around the true value for small $n$, but the discrepancy decreases as $n$ increases.
This simple exercise provides a gentle introduction to Monte Carlo methods: simulation allows us to approximate theoretical quantities, and the approximation improves systematically with larger sample sizes.
\paragraph{Standard error.} Let $I_i = \mathbf{1}\{X_i < -2\}$ where $X_i \sim N(0,1)$. Then $I_i \sim \text{Bernoulli}(p)$ with $p = \Pr[X < -2]$. The empirical probability is the sample mean \[ \hat{p} = \frac{1}{n}\sum_{i=1}^n I_i. \] Since $\mbox{Var}(I_i) = p(1-p)$, the variance of $\hat{p}$ is $mbox{Var}(\hat{p}) = p(1-p)/n$. Hence the standard error is \[ \text{SE}(\hat{p}) = \sqrt{\tfrac{p(1-p)}{n}}. \] In practice $p$ is unknown, so we replace it by $\hat{p}$: \[ \widehat{\text{SE}} = \sqrt{\tfrac{\hat{p}(1-\hat{p})}{n}}. \]
Consider the task of approximating the value of $\pi$. Imagine randomly throwing darts at a square board with a circular target inscribed inside. If you know the total number of darts thrown and the number that land inside the circle, you can use these quantities to estimate the value of $\pi$. This is because the ratio of the area of the circle to the square approximates $\pi / 4$. By throwing a large number of darts (or, in a computational sense, generating a large number of random points) and calculating the proportion that fall inside the circle, you can multiply by $4$ to get an approximation of $\pi$.
We define a circle of radius $r=1$. Then its surface is known to be $\pi r^2$. We draw $x$ and $y$ from the uniform distributions $a=-1$ and $b=1$. We write the probability density function of the uniform distribution:
The square surface will be $2\times 2 = 4$ and the area enclosed by the circle is $\pi$:
A recurring question is on the choice of the number of draws for the Monte Carlo method. We are not aware of a definite answer, but here is an illustration of the histogram of $100$ estimates of $\pi$ we found with $10^4$ draws:
compared with $100$ $\pi$ estimates we found with $10^6$ draws:
\stepcounter{question}
\paragraph{Remark (Monte Carlo integration).}\footnote{For an accessible introduction, see \href{https://www.scratchapixel.com/lessons/mathematics-physics-for-computer-graphics/monte-carlo-methods-in-practice/monte-carlo-integration}{Scratchapixel, Monte Carlo Integration}.} Monte Carlo integration interprets the value of an integral as an expected value. If $X$ is uniformly distributed over a domain $\mathcal{D}$ with volume $V$, then \[ \int_{\mathcal{D}} f(x)\,dx \;=\; V\,\mathbb{E}[f(X)]. \] The corresponding estimator with $N$ independent samples $\{X_i\}_{i=1}^N$ is
It is unbiased, with variance \[ \operatorname{Var}(\hat{I}_N) \;=\; \frac{V^2}{N}\,\operatorname{Var}(f(X)). \] Thus the error decreases at rate $N^{-1/2}$, regardless of the dimension of $\mathcal{D}$.
\stepcounter{question}
\paragraph{Example (estimating $\pi$).} In Section (ref) we approximated $\pi$ using the area of a quarter circle. Define the indicator function \[ H(x,y) =
\] Sampling uniformly on the square $\mathcal{D}=[-1,1]\times[-1,1]$ of volume $V=4$, the probability that $(x,y)$ falls inside the unit disk is $\pi/4$. The Monte Carlo estimator is therefore \[ \hat{\pi}_N \;=\; V\,\frac{1}{N}\sum_{i=1}^{N} H(x_i,y_i). \] As $N\to\infty$, $\hat{\pi}_N$ converges to $\pi$ with variance decreasing at rate $1/N$.
\paragraph{Estimate $\displaystyle \int_0^1 e^{x}\,dx$.} Let \[ I \;=\; \int_{0}^{1} e^{x}\,dx \;=\; e-1. \] The Monte Carlo estimator based on i.i.d.\ draws $U_1,\dots,U_N \sim \mathcal U[0,1]$ is
\paragraph{Unbiasedness and LLN.} Since $\mathbb{E}[e^{U_i}] = \int_0^1 e^{x}\,dx = I$, we have $\mathbb{E}[\widehat I_N]=I$. By the law of large numbers, \[ \widehat I_N \xrightarrow[]{\;\text{a.s.}\;} I \quad \text{as } N\to\infty. \]
\paragraph{Variance and CLT.} With $U\sim \mathcal U[0,1]$, \[ \operatorname{Var}(\widehat I_N) \;=\; \frac{1}{N}\,\operatorname{Var}\!\big(e^{U}\big). \] Compute \[ \mathbb{E}\!\left[e^{tU}\right] \;=\; \int_0^1 e^{tx}\,dx \;=\; \frac{e^{t}-1}{t}, \] so \[ \mathbb{E}\!\left[e^{U}\right]=e-1, \qquad \mathbb{E}\!\left[e^{2U}\right]=\frac{e^{2}-1}{2}, \] and therefore \[ \operatorname{Var}\!\big(e^{U}\big) \;=\; \frac{e^{2}-1}{2}-(e-1)^2 \;\approx\; 0.2420356. \] Hence \[ \operatorname{Var}(\widehat I_N) \;=\; \frac{0.2420356}{N}, \qquad \operatorname{se}(\widehat I_N) \;=\; \sqrt{\frac{0.2420356}{N}} \;=\; \mathcal O\!\left(N^{-1/2}\right). \] By the central limit theorem, \[ \sqrt{N}\,\big(\widehat I_N - I\big)\;\xrightarrow{d}\; \mathcal N\!\left(0,\operatorname{Var}\!\big(e^{U}\big)\right), \] so an approximate $(1-\alpha)$ confidence interval is \[ \widehat I_N \;\pm\; z_{1-\alpha/2}\,\sqrt{\frac{\widehat{\operatorname{Var}}(e^{U})}{N}}, \] where $\widehat{\operatorname{Var}}(e^{U})$ is the sample variance of $\{e^{U_i}\}_{i=1}^N$.
One of the most important applications of Monte Carlo methods is the approximation of high-dimensional integrals. Such problems arise in Bayesian statistics, statistical physics, and quantitative finance, where analytic solutions are not available.
\paragraph{Bayesian example.} Consider data $D$ and a parameter vector $\theta \in \mathbb{R}^{p}$. The posterior distribution is \[ \pi(\theta \mid D) \;=\; \frac{L(D \mid \theta)\,\pi(\theta)}{\int_{\mathbb{R}^{p}} L(D \mid \vartheta)\,\pi(\vartheta)\,d\vartheta}, \] where $L(D \mid \theta)$ is the likelihood and $\pi(\theta)$ the prior. The denominator,
is called the marginal likelihood or model evidence. For $p \geq 5$, this integral has no closed form and is effectively impossible to evaluate with deterministic numerical quadrature, since the cost grows exponentially with $p$ (the “curse of dimensionality”).
\paragraph{Monte Carlo solution.} Using samples $\theta^{(1)},\dots,\theta^{(N)}$ drawn from the prior $\pi(\theta)$, one may approximate (ref) by \[ \hat{Z}_N \;=\; \frac{1}{N}\sum_{i=1}^N L(D \mid \theta^{(i)}). \] This is an unbiased estimator of $Z$, with variance decreasing at rate $1/N$. More sophisticated methods (importance sampling, Markov chain Monte Carlo) improve efficiency, but all rely on the same core principle: Monte Carlo integration is often the only practical tool to evaluate otherwise intractable high-dimensional integrals.
Monte Carlo methods are not only useful in statistics but also in finance. Many derivative securities are priced by computing expectations under a risk–neutral measure, often involving high-dimensional stochastic processes with no closed-form solution.
\paragraph{Inflation cap.} An inflation cap is a contract that pays off if the realized inflation rate $\pi(n)$ over $n$ years exceeds a fixed strike $k$. The payoff at maturity is \[ \text{Caplet Payoff} \;=\; \max\!\big(\pi(n)-k,\,0\big). \] The present value at time 0 under the risk–neutral measure $\mathbb{Q}$ with discount factor $B(0,n)$ is
\paragraph{Intractability.} In practice, $\pi(n)$ is modeled as a stochastic process (e.g. with stochastic volatility or jumps). The distribution of $\pi(n)$ has no closed form, so the expectation in (ref) cannot be computed analytically.
\paragraph{Monte Carlo solution.} We simulate $N$ independent paths of the inflation process under $\mathbb{Q}$, producing realizations $\pi^{(1)}(n),\dots,\pi^{(N)}(n)$. The Monte Carlo estimator of the caplet price is \[ \widehat{C}_N \;=\; B(0,n)\,\frac{1}{N}\sum_{i=1}^N \max\!\big(\pi^{(i)}(n)-k,\,0\big). \] This provides a consistent and unbiased approximation of (ref), with accuracy improving at the standard Monte Carlo rate $1/\sqrt{N}$.
\paragraph{Modeling inflation dynamics.} To simulate $\pi(n)$, we need a model for the inflation process under the risk--neutral measure $\mathbb{Q}$. In practice, one proceeds as follows:
We suggest exercises\footnote{pandas_exercise_part1.py} to get familiar with some pandas basics features in python's pandas library.
First we want to import the package panda and create a dataframe we call df:
You might want to sum over rows for each column .sum(axis=0) or sum over columns for each row .sum(axis=1)
The panda's documentation explains how to access specific element of a DataFrame, using either the label or the integer positions. \href{http://pandas.pydata.org/pandas-docs/stable/indexing.html#selection-by-label}{Link here}.
You can use the .loc attribute to access elements in a DataFrame, label-based. For df, that would be:
You can also create subset of your dataframe, slicing it.
Which yields from left to right:
If you want to be able to select a cell in a DataFrame, regardless of the name of the column or the name of the row, you can use integer based indexing with .iloc attribute. Keep in mind that python indexing convention starts at $0$.
Creating the above slices would become:
Note the subtlety, illustrated by our example slice3i:
You can also use other selection methods, described in panda's documentation: .at, .iat, .ix
We can use the function unstack on the DataFrame to transform it into a wider format where multiple levels of data are reshaped based on combinations of row and column indices. This is especially handy for panel regressions, as it aids in disentangling data across multiple entities (like individuals or firms) and time periods, providing a clearer view of observations over time. For the given DataFrame, when we apply unstack, we observe for df_unstacked:
By using positions or names of columns and/or index, we are able to loop through a DataFrame.
The idea behind this is that the elements of the DataFrame df2 are found:
\quad
\quad
It is also possible to transpose a data frame. The transpose of df2 would be:
\quad
We want to create\footnote{pandas_exercise_part1_ts.py} the following DataFrame with the French population values:
\quad
Look carefully at the way your DataFrame is indexed, it starts from 0 and not from 1. Therefore here you have 9 observations indexed from 0 to 8.
As we entered the data in the wrong order, we need to reverse:
You can look into python's documentation on slice class. Imagine you have a list you call mylist:
Conveniently, we chose the value of each item to be its index position in the list, so $0$ is at index $0$, $1$ at index position $1$, etc. Imagine you want to reverse this mylist: we have a list of length $4$ so going from index position $0$ to index position $3$. To enumerate the index, you can do: $start:end:step$ within the $[$ and $]$ signs. So mylist[0:len(mylist):1] will return the full list.
To reverse the list, we want the step to be $-1$, we omit the start and end position, so that python takes all possible items. Just look back at section (ref) where we used $1:$ to select all items starting at index $1$ and after.
This is optional, but you might want to index your DataFrame with time. First you want to create an object with the dates as 'year-month' from 2016-01 to 2016-09 and modify our DataFrame index.
Note that we have to range until 2016-10 for the last item of our dates object to be 2016-09.
We see that the name of our column is '0' we want to rename it with 'population'. We can have a look at the data by plotting it and choosing a title.
We want to compute the change of the French population from month $m$ to month $m+1$: $$\mbox{change}_{m+1}=\frac{\mbox{pop}_{m+1}-\mbox{pop}_m}{\mbox{pop}_m}$$
In visual terms, that will be:
\quad
Arthur Cayley worked on matrices and introduced operations cayley1858ii.
We create\footnote{pandas_exercise_part1_matrices.py} a vector with its index as the matrix column names and we can then perform a matrix multiplication\footnote{It is important that the columns names of the matrix are the same as the index of the vector for the product to be acceptable}
\[
\times
=
\]
We might want to multiply element by element, this is also called Hadamard product, that we write here $\odot$.
\quad
First we duplicate the column to form a square matrix and multiply the two matrices:
\[
\odot
=
\]
A matrix is fundamental to many advanced concepts in mathematics and econometrics.
Another way to determine if a matrix is non-singular is to look at its rows (or columns). If all rows (or columns) of an $n\times n$ matrix are linearly independent -- meaning no row (or column) can be expressed as a combination of the others -- then the matrix is non-singular. This is also referred to as the matrix having "full rank". A quick check for this is the determinant: non-singular matrices have a determinant that is not zero.
In the realm of econometrics, matrices play a vital role in various models, including AR(p) and VAR models.
When analyzing the stability of a linear propagation process, it is useful to study the eigenvalues and eigenvectors of the underlying matrix.
An eigenvector $v$ of a matrix $\Delta$ is a non-zero vector such that \[ \Delta v = \lambda v, \] where $\lambda$ is the associated eigenvalue. Intuitively, multiplying $\Delta$ by $v$ only stretches or shrinks $v$, without changing its direction.
This matters because the long-run behavior of repeated applications of $\Delta$ depends on its eigenvalues. If all eigenvalues are strictly smaller than one in absolute value, repeated multiplications of $\Delta$ dampen shocks. If at least one eigenvalue exceeds one, shocks are amplified and the system becomes unstable.
\paragraph{Application to financial networks: DebtRank} The matrix $\Delta$ summarizes how distress propagates from one bank to another. Bardoscia2015aa model contagion by transmitting only the newly incurred distress at each step. The update rule is \[ h(t+1) = h(t) + \Delta \big(h(t)-h(t-1)\big). \] Here $h_i(t)$ represents the relative equity loss of bank $i$ at time $t$. When bank $j$ suffers an additional equity loss, the counterparties $i$ that are exposed to $j$ must reduce the value of their interbank assets accordingly. This is captured by the matrix $\Delta$, which scales the shock according to the size of exposures relative to each bank's equity base. In this way, distress is linked directly to a deterioration of capital positions and thus to the valuation of exposures. Each bank's distress is capped at $h_i(t)\leq 1$, reflecting the fact that no institution can lose more than 100% of its equity.
Suppose each bank is hit at time $t=1$ by an initial shock $h(1)=0.5\%$. At this stage, the losses are simply $h(1)$.
At $t=2$, the first round of contagion takes place: the initial shock is transmitted once through the network, giving an additional contribution $\Delta h(1)$. The cumulative distress is then \[ h(2) = h(1) + \Delta h(1). \]
At $t=3$, the second round of contagion occurs: the shock passes along paths of length two, adding $\Delta^2 h(1)$. The cumulative distress is \[ h(3) = h(1) + \Delta h(1) + \Delta^2 h(1). \]
Continuing in this way, after $t$ steps we have \[ h(t) = \sum_{s=0}^{t-1} \Delta^s h(1). \]
We have \[ S_t = I + \Delta + \Delta^2 + \cdots + \Delta^{t-1}. \]
For scalars you know that \[ 1 + r + r^2 + \cdots + r^{t-1} = \frac{1-r^t}{1-r}, \quad |r|<1. \]
The matrix version is analogous: \[ S_t = (I - \Delta^t)(I-\Delta)^{-1}, \] provided that $I-\Delta$ is invertible.
If all eigenvalues of $\Delta$ are strictly smaller than one in modulus, then the infinite sum converges and the long-run impact of the shock is finite: \[ h(\infty) = (I-\Delta)^{-1} h(1), \] with $I$ the identity matrix. If instead at least one eigenvalue exceeds one, the sum diverges and the shock is amplified without bound.
Cholesky decomposition is a powerful mathematical tool often used in econometrics.
Cholesky decomposition:
For a matrix $A$ that is a real positive-definite, the Cholesky decomposition breaks it down into two components. It represents $A$ as the product of a lower triangular matrix, $L$, and its transpose, $L^T$. So, the decomposition can be represented as: $A = L L^T$
What's striking about this is that the matrix $L$ is unique, with real and positive diagonal entries. This uniqueness is paramount, especially when dealing with SVARs in econometrics, as it ensures consistent and unambiguous results.
As a reminder of matrix and random variables: if $X$ is a collection of rvs such that its mean is $\mu$ and its variance-covariance matrix is $Var(X) = \Sigma$. Now if we take a matrix composed of constant elements $U$ and define a new random variable vector $Y=UX$, then $E(Y)=U \mu$ and $Var(Y)=U Var(X) U^T = U \Sigma U^T$.
A Cholesky decomposition can be used useful to generate multivariate normally distributed random vectors. If $Z=[Z_1,\ldots,Z_N]^T$ is a vector where $\forall i, Z_i \sim N(0,1)$ and are independent, then if $U$ is a $(N\times N)$ matrix,
If $C$ is the covariance matrix we are looking to introduce, then it is a real symmetric positive-definite matrix and its Cholesky decomposition is $C = U^T U$.
If we follow the spirit of Kousky:2012wg and the lecture of \href{https://martin-haugh.github.io/files/MonteCarlo/MCS_Generate_RVars.pdf}{Martin Haugh}. I generate $3$ iid rvs $Z$ from a standard normal law, with samples of size $10^5$. Then I want to make those rvs dependent with a covariance matrix: \[ C =
\]
Then I extract $U^T$ from the Cholesky decomposition of $C$ and multiply $U^T Z$
Computational complexity\footnote{pandas_exercise_part1_cc.py} concerns the resources required to solve computational problems, primarily time and space. The theory classifies algorithms according to their performance in relation to the input size. Let us denote the input size by $n$.
To begin, let us define the Big O notation. For a given function $g(n)$, we say that $f(n) = O(g(n))$ if there are positive constants $c$ and $n_0$ such that $0 \leq f(n) \leq c \cdot g(n)$ for all $n \geq n_0$. This notation helps us describe an upper limit on the growth rate of the runtime or space requirements of an algorithm.
Consider a simple example: an algorithm that checks each pair in an array to find duplicates. This algorithm would have a time complexity of $O(n^2)$ because, in the worst case, it performs a comparison for each pair of elements in the array.
P (Polynomial Time): This class contains problems that can be solved by a deterministic Turing machine in polynomial time. An example is matrix multiplication, which can be solved in time $O(n^3)$ using the standard algorithm, though faster algorithms exist.
Consider two $2 \times 2$ matrices $A$ and $B$ defined as follows: \[ A =
, \quad B =
\] The product of these matrices, $C = A \cdot B$, is calculated as: \[ C =
\] Each element of matrix $C$ is the sum of products of elements from a row of $A$ and a column of $B$. For a general $n \times n$ matrix, this calculation involves $n$ multiplications and $n-1$ additions for each of the $n^2$ entries, resulting in a total of $n^3$ operations, assuming the addition and multiplication have a constant time complexity. Hence, the naive matrix multiplication complexity is $O(n^3)$.
The Strassen algorithm, introduced in strassen1969gaussian, reduces the computational complexity of matrix multiplication. This algorithm is particularly effective for large matrices and is based on a divide-and-conquer approach that reduces the number of recursive multiplications.
Consider the same $2 \times 2$ matrices $A$ and $B$, instead of performing eight multiplications as in the naive approach, the Strassen algorithm performs only seven multiplications by introducing intermediate products:
These intermediate products are used to compute the entries of the product matrix $C$ as follows: \[ C =
\]
By reducing the number of multiplications, the Strassen algorithm improves the overall time complexity. For general \(n \times n\) matrices, the algorithm recursively applies itself to smaller matrices, leading to a complexity of \(O(n^{\log_2 7}) \approx O(n^{2.807})\). This improvement is significant, particularly for very large matrices.
This topic is still improved, as in duan2023faster and williams2024new.
We suggest\footnote{pandas_exercise_part1_import.py} to build a DataFrame of the evolution of the French population as taken from INSEE website: \href{http://www.insee.fr/fr/bases-de-donnees/bsweb/serie.asp?idbank=001641607}{here}.
Now you need to 'tell' python where you file is located.
Import the package Miscellaneous operating system interfaces:
Then you need to tell him where your data folder is located, in my computer, that would be:
If you are note sure where your file is located, here is a step-by-step guide:
\quad
Another way to approach this is by using the Spyder working directory browser until you find the folder where your file is located:
From the INSEE website \href{http://www.insee.fr/fr/bases-de-donnees/bsweb/serie.asp?idbank=001641607}{here} you can choose to import the table as a csv file. It will be called 'valeurs.zip'.
We prefer to work with comma-separated values file, for some of the reasons described \href{https://en.wikipedia.org/wiki/Comma-separated_values}{here}.
First unzip the file to have it a Valeurs.csv
Copy the Valeurs.csv file into your usual data folder.
Follow the indications section (ref) to indicate where your file is located.
If you open your file with a text editor\footnote{be carefule, if you open this file with Excel, do not save the modifications}, it looks like this:
So because we are working with French data:
Hence we need to tell that:
The data set that has been imported and is stored into the DataFrame named df is not clean.
First we want to have an timely ascending series:
We want to rename the columns of the DataFrame with Month and Population.
We want to reset the index, knowing that the observations are monthly and start in January 1994:
Now that we have reset the index, we can drop the column with the months that is not bringing much information now, you drop vertically, a column, that is axis $= 1$.
Now try to plot your DataFrame df, it doesn't work and tells you:
This means that when importing the data from the .csv file, python did not recognise that the column 'Population' were numeric. We now need to convert those data to numeric. First you need to remove the space in the numbers, then converting strings to float (note the use of regex\footnote{Without regex=True the method will look for exact matches. When you use regex=True it will look for the sub strings too}):
Now when you try to plot, you should get:
There is the possibility to resample our DataFrame from one frequency to another using the resample() function. We need to indicate the desired frequency described in the Offset Aliases section of the time series documentation available \href{http://pandas.pydata.org/pandas-docs/stable/timeseries.html}{here}. We list some of the more commonly used:
The we need to chose a "method" for resampling, for example .sum() or .mean()
In our case, if we want to move to yearly frequency and get the average over the 12 months. If we plot the obtained DataFrame, we see that the plot is slightly smoother.
Now you can plot and compute monthly changes.
You can get some descriptive statistics on the population and the monthly population changes:
Wassily Leontief is widely known for developing the Input-Output model of the economy. His foundational work on the subject began with Leontief:1949uo. For a comprehensive collection of his major papers on this model, one can refer to the compilation presented in the book leontief1986.
The essence of the model lies in its sectoral division of the economy. Every sector is dependent on others for its input to produce its output. Thus, the production of a sector isn't just shaped by external demand (from end consumers) but also by internal demand (from other sectors).
To crystallize this concept using matrix notation, let's define the following:
With these definitions in place, the relationship between internal and external demands in determining sectoral production is given by:
Rearranging the above equation, we get:
Here, \(I\) represents the identity matrix. If the matrix \(\left(I-M \right)\) is invertible, we can confidently assert that our equation has a unique solution equation (ref). On the contrary, if it's non-invertible, the model becomes indeterminate, meaning there could be no solution or potentially infinite solutions.
Let's imagine we have the simplified economy:
We use the consumption matrix $M$ after dividing by the Total gross output. We can now estimate following equation (ref)
Now, if there is a one unit increase in the demand for agriculture output, this should have a positive impact of $.2$ on the manufacturing production.
This is a step-by-step way of doing it\footnote{input_output_illustration.py} Actually, the most difficult part is to divide each column of a matrix by a given value from a vector, three approaches are suggested (as seen during class):
\stepcounter{question}
\stepcounter{question}
\stepcounter{question}
I would tend to say that before engaging in a long research project, you should explore whether some data will be available to test your model or not.
You will need to change the working directory to be able to load some files,previous section\footnote{get_data_vansteenberghe.py} explained step-by-step how to change your working directory.
Banque de France provides historical times series via its webstat portal, \href{http://webstat.banque-france.fr/fr/}{here}.
We might want to download the historical evolution of property prices in France, this can be found (with some more explanations) under the code RPP.Q.FR.N.ED.00.1.00 (Indices des prix des logements anciens, Ensemble des logements, France metropolitaine).
\stepcounter{question}
You can access INSEE data from \href{https://www.insee.fr/fr/statistiques}{their website}. As an illustration, the Demography - Population at begining of a month in France, series ID: 000436387
\stepcounter{question}
The OECD provides panel data via its websilte \href{https://stats.oecd.org/}{here}.
You can select "index of Industrial Production". Then you'll need to reorganise the data set, making use of pivot table:
You can get data from the Federal Reserve of Saint Louis, their Federal Reserve Economics Database, or FRED.
What is convenient is that you can download the time series directly form python if you know the code of the variable you are looking for. For example, we can download the following variables for the US economy, which are the usual variable used in basic macroeconomic models:
\stepcounter{question}
Sometimes, it is possible to get data directly from providers like Euronext, where it is possible to get index historics \href{https://www.euronext.com/en/indices/index-statistics?archive=this_year}{here}. This enables you for example the return of the CAC index with or without the dividends:
With \href{https://www.quandl.com/tools/python}{Quandl, you can "get millions of financial and economic datasets from hundreds of publishers directly into Python"}. You can start by \href{http://help.quandl.com/article/93-how-do-i-install-the-quandl-python-library}{installing} the Quandl library\footnote{if you are using Windows, you might have to use a Command Prompt: C:/Anaconda3/Scripts/pip install quandl}, then you can \href{http://help.quandl.com/article/205-how-do-i-download-a-dataset-using-python}{download datasets}.
You should request your API by creating an account. You can then download massive dataset, \href{http://help.quandl.com/article/92-how-do-i-download-the-quandl-codes-of-all-the-datasets-in-a-given-database}{as explained here}.
In order to get financial time series from Yahoo finance, you first need to define the stock (or list of stocks) you are interested in. For this, you can search at the Yahoo finance URL, \href{https://finance.yahoo.com/}{here} for the name of the underlying company, let's say for illustration purpose Deutsche Bank.
You get historical time series on opening and closing prices, highs and lows as well as volumes.
At the University, you have access to a Bloomberg terminal. We suggest that you install the Bloomberg extension for Excel (look into the start menu in window, you should find this installer).
\quad
You might want to download some historical time series with Excel. We suggest the following step:
You need to go onto the Bloomberg terminal in order to find the correct field: type FLDS
Then you can search for a company name for example, it will give you its Bloomberg ticker:
Once you get the Bloomberg ticker (e.g. "KN FP Equity" for the bank BPCE), you can use this ticker in an excel spreadsheet.
We suggest to organize your Excel sheet as follow:
Here is an simple example to download stock price data vertically:
You can get data from the IMF \href{http://data.imf.org/?sk=388DFA60-1D26-4ADE-B505-A05A558D9A42}{here}, from the ECB at its \href{http://sdw.ecb.europa.eu/}{Statistical Data Warehouse}. The World Bank provides time series on \href{https://data.worldbank.org/data-catalog/world-development-indicators}{World Development Indicators}. The Bank for International Settlements (BIS).
\stepcounter{question}
We do not intend to provide here a full lecture on web scraping with python. We provide one use case: getting a list of Fintech that are IFP from \href{https://www.orias.fr/web/guest/search}{ORIAS website}.
Lending in France is regulated by law under what is known as monopole bancaire\footnote{Article L511-5, this monopoly is a mean to protect depositors from institutions insolvency and illiquidity} and is regulated at the European level by the CRDIV\footnote{Article 9-2 of the directive 2013/36/EU of the European Parliament and of the Council of 26 June 2013 on access to the activity of credit institutions} on taking deposits or other repayable funds and by the CRR\footnote{Article 4-1 (1) of the regulation (EU) No 575/2013 of the European Parliament and of the Council of 26 June 2013} for granting credits for its own account. This regulation is a barrier to entry into the market but has been relaxed by the introduction of new categories of actors in France in 2014\footnote{Ordonnance n° 2014-559 of the 30th May 2014 which is a Dérogation au monopole bancaire introducing the role of conseillers en investissements participatifs and intermédiaires en financement participative}. For the French market, the ACPR\footnote{Autorité de contrôle prudentiel et de résolution, the French financial regulator} has the authority to deliver accreditation\footnote{Ordonnance n° 2014-559 of the 30th May 2014 article 16} to intermediaries in participatory financing (IFP). ORIAS provides information on the IFP but to gather the information, one has to research and select the data. We suggest to use python libraries for web scraping to feed a data frame:
\resizebox{\textwidth}{!}{
}
This exercise was done with the help of a StackOverflow member (for more details on the problematic see \href{https://stackoverflow.com/questions/57536576/python-3-how-to-scrape-research-results-from-a-website-using-csfr/57539377#57539377}{StackOverflow website}) that helped us overcome the CSFR restrictions.
We provide the code: 20190817_orias_scrape.py with little comments, this is for advanced python users, you might skip this part for now and come back on it later on.
\stepcounter{question}
In the domain of Machine Learning, as well as in econometrics, data assurance quality—encompassing activities like data cleaning, preparation, and organization—typically constitutes a staggering $80\%$ of the overall workload. Beyond just the effort involved, ensuring data quality is pivotal. It not only enhances the accuracy and reliability of the model but also furnishes insights into the inherent strengths and potential vulnerabilities of the dataset.
ottonello2020financial works on a panel of firms:
To elucidate this crucial aspect, we\footnote{vansteenberghe_data_quality_ML_US_wages.py} turn our focus to the US Adult Census dataset, which correlates income with sociodemographic attributes such as age, education, race, and other pertinent features. This dataset is grounded in the seminal work of kohavi. For those seeking to delve deeper, the \href{https://archive.ics.uci.edu/ml/datasets/census+income}{data set can be accessed here}.
An indispensable starting point is to familiarize oneself with the dataset and promptly identify any missing values:
A pertinent question arises: are the NaNs indicative of non-working individuals, or are they the result of data collection errors? If it's the latter, the optimal course of action might be to discard the affected rows.
Upon examining the 'native_country' feature, an anomaly surfaces—a country labeled "South". A detailed inspection, especially when cross-referenced with the 'race' attribute, strongly suggests that this likely corresponds predominantly to South Korea.
Feature engineering can serve as a potent tool to enhance the data's utility. For instance:
We\footnote{vansteenberghe_ML_wage.py} did some preliminary work as a starting point. \stepcounter{question}
Section (ref) we explored in great details a parametric model $m$, linear and at best quadratic for a variable of interest ($y$: log of wage) using explanatory variables also called features ($X$: gender, education, experience) in $y = m(X) + \epsilon$. We might wonder why not to use machine learning for our regression problem.
We suggest an answer following charpentier2018, considering that econometrics (statistics)
, to make inference, we are making assumptions about the data generating process (cf. section (ref)). If we had to distinguish machine learning from econometrics, we would say that machine learning is focused on predictive power and does not make any assumption on the data generating process.
In both econometrics and machine learning, having data sets $[y_i,X_i]$ the goal is to "search" for a model that minimizes a loss function:
under some constraint on $m$.
We\footnote{The main code to follow is ML_credit_cards.py} present some machine learning technics\footnote{This is inspired by the lectures of Emmanuel Rachelson \href{http://emmanuel.rachelson.free.fr/}{(his website)} and the book raschka} applied to some data set with the python code.
In the words of Andrews:2017vv:
and their work expand this to nonlinear econometric models. Such transparency is not straightforward for machine leaning models.
\href{https://www.youtube.com/watch?v=jBxfQJOHJxc}{Some French reporter investigated the testing methods used by the Union Cycliste Internationale (UCI)} to test whether there were hidden motors in bicycles before a race. The reporter questioned the efficacy of the methods used... Could it be that there are engines used during the Tour de France?
The UCI \href{https://www.sportsintegrityinitiative.com/uci-statement-technological-fraud/}{replied with a published note}.
What is very interesting, is that reading this written reply, I started to have strong doubt... The UCI methods could be inefficient.
The UCI is quoting a report performed by a laboratory in the United States of America, the Microbac report, \href{https://www.uci.org/docs/default-source/imported-library/uci17d0410_rp1_neutral.pdf?Status=Master&sfvrsn=f0e6cf5b_6}{available here, under the code name UCI17D0410_rp1_Neutral.pdf}.
What is quite disturbing is that the report is actually what led me to doubt the methods at UCI:
When the detection tool is used on a bike, at 10 millimeters from the bike frame, there are 7 instance of "10" out of $4 \times 4 = 16$ tests. A "10" means that you have a strong signal, a positive. Out of those positive, 4 where bike with actual engines (in yellow)... but 3 of those "10" had no engine inside.
\quad
In other terms, testing the bikes, the Microbac found 3 false positive out of 16 tests, that is a ratio of $\frac{3}{16}$.
\quad
If you read carefully the UCI answer, they performed "over 4,000 controls undertaken during the 2017 Tour de France", this means that --if the Microbac tests have any meaning-- the UCI testers have found an expected 750 positive signals (all false positive, as no scandal arose from a true positive engine in a bike during the Tour de France 2017).
\quad
The question I would have is: what did the UCI do of those close to 1,000 positive alarms during the 2017 Tour de France? Did they take the time to perform thorough testing (computer tomography)... or just let the biker race?\footnote{I actually sent an e-mail to the UCI team, they answered along the lines that using the test device isn't just as simple as seeing a 10 on the screen and assuming it is a positive. The size and shape of the reading on the test device screen also needs to be interpreted which is why they deliver 3-4 hours of training to each operator. They claimed that they would not get many false positives because there is more to using the device than just seeing the number on the screen. They would occasionally get a positive reading that appears to be false and dismantle the bikes but they replied that it is very rare and not statistically significant.}
A model becomes inefficient if it yields a disproportionately high ratio of False Positives. In such scenarios, there might be a need to manually differentiate between True and False Positives.
For supervisory entities:
A relevant case to consider is the Dieselgate scandal. This matter is of high importance as it pertains directly to public health. Were testing protocols for this consistently the same? Did it never occur to testers to sporadically test vehicles on roads? Was their incentive structure built around identifying dishonest auto manufacturers? Did organizations like the UCI genuinely have an interest in identifying deceitful racers, even if it meant another public scandal?
These incentive and testing protocol inquiries are equally relevant for financial supervisors. Additionally, testing protocols should be objective and not constructed based on a desired result.
One certainty is that when using train services, the companies appear motivated to identify and penalize fare evaders. They maintain secrecy about their checking methods, including the frequency, methodology, and locations of inspections.
We will first train our model on the Train data set, then use the model to predict the default $Y$ on the Test data set. We want that our model predict correctly default and non-default: $Y_{pred}=Y$, we define true positive, true negative, false positive and false negative as follow:
\quad
From the confusion matrix, we can define several metrics, here we focus on three:
False Positive Rate: $FRP=\frac{FP}{Total}$
Classification Accuracy: Overall, how often is the classifier correct? $ CA = \frac{TN + TP}{Total}$
Sensitivity: When the actual value is positive, how often is the prediction correct? $S = \frac{TP}{FN + TP}$
With our model, we are attributing two probabilities, one for each state : default ($0$) and non-default ($1$).
For illustration, if with our model we have Proba($0$) > Proba($1$) but $Y=1$, then our model doesn't predict a default while in fact a default was observed, it is a false negative.
The Receiver Operating Characteristic (ROC) curve hanley1982meaning,fawcett2006introduction is a graphical representation used to evaluate the performance of a binary classifier system as its discrimination threshold is varied. Specifically, the ROC curve plots the True Positive Rate (sensitivity) against the False Positive Rate (1-specificity) for different threshold values. It serves as a tool to illustrate the trade-offs between sensitivity and specificity for a given classifier.
We also define:
The scikit learn community provides a comprehensive guide to machine learning techniques. The first question one would have when faced with a problem is the way to choose the estimator. \href{http://scikit-learn.org/stable/tutorial/machine_learning_map/index.html}{They provide a flow chart on their website}.
Scikit-learn also provide a \href{http://scikit-learn.org/stable/user_guide.html}{user guide that is quite exhaustive}.
We will mainly work on the data set introduced section (ref) and as in the paper Yeh2009aa, we want to test the following data mining techniques:
We add to that list:
We\footnote{ML_wage_naive_bayes.py} use the data set from Schooling, Experience and Earnings of Mincer 1974 work. We find a \href{https://archive.ics.uci.edu/ml/machine-learning-databases/00350/default
\quad
\quad
To determine the sex of observation 6, let's denote the sex as \(S\), where \(S=1\) stands for male and \(S=0\) for female. The variables \(W\), \(A\), and \(E\) represent Wage, Age, and Education, respectively.
Our primary objective is to determine if observation 6 is male. Formally, this can be expressed as: \[ \mathcal{P}\left(S=1 \mid (W,A,E)\right) > \mathcal{P}\left(S=0 \mid (W,A,E)\right) \]
To determine the above probabilities, we can use Bayes' theorem, which states:
Here:
From our sample, we know that: \[ \mathcal{P}\left(S \right) = 54\% \]
If we discretize Wage, Age, and Education into \(x\) segments, we then create a 3-dimensional array which results in \(10^x\) values. This expansion of values due to multiple dimensions is often termed as the 'curse of dimensionality'.
The conditional probability of observing \(W\), \(A\), and \(E\) given the sex can be represented as: \[ \mathcal{P}\left( (W,A,E) \mid S\right) = \mathcal{P}\left( W\mid S\right) \times \mathcal{P}\left( A \mid (S,W)\right) \times \mathcal{P}\left( E \mid (S,W,A)\right) \]
However, to simplify calculations and reduce the complexity, we can make the naive assumption that wage, age, and education are conditionally independent given the sex. Therefore: \[ \mathcal{P}\left( (W,A,E) \mid S\right) = \mathcal{P}\left( W\mid S\right) \times \mathcal{P}\left( A \mid S\right) \times \mathcal{P}\left( E \mid S\right) \]
With this assumption, we are left to estimate only \(3 \times x\) values.
Nota bene: In equation (ref), we don't need to compute the denominator. It serves as a normalizing constraint to ensure that the probabilities of being male or female sum to one, given that we consider only binary categories for individuals.
We make the assumptions that Wage, Age, and Education are normally distributed, once we have the mean and the standard deviation, then the probability can be simply derived using: $$\mathcal{P}\left( x\vert (\mu , \sigma )\right) = \frac{1}{\sqrt{2 \pi \sigma^2}} e^{-\frac{(x-\mu)^2}{2 \sigma^2}}$$
We obtain the following confusion matrix:
Which yields a F1 score of $0.65$, we tend to accept this rather simple model.
KNN is a "lazy" learner, meaning that it doesn't learn a discriminative function but rather memorize the training data set.
The KNN algorithm follows the following logist and steps:
In this illustrative example, $k=6$ and we use the Euclidean distance (assuming a two dimension problem). The value of $x$ will be $2$, by majority.
If we\footnote{ML_wage_KNN.py} use again Mincer's 1974 work.
We take a random sample of $21$ elements from the data set. We change the last observation SEX value to $99$ and we want to predict its value ($0$ or $1$) using a 5-nearest neighbors method.
We then compute the Euclidean distance between this point ($99$) and the rest of that sample. We then count the SEX of the $5$ nearest neighbors and decide on the SEX of our point based on majority vote.
In a two dimensions example the Euclidean distance is simply: $\sqrt{(x_2-x_1)^2 + (y_2-y_1)^2}$.
If we\footnote{ML_wage_DTC.py} use again Mincer's 1974 work.
The obtained Decision Tree classifier is:
In our case, we have limited the depth of our tree to $3$. This divide our 2-Dimension space as follow:
In practice, we are using a binary tree, where each parent has two children. We can in this case simply define the entropy as: $$Ent= - p(l) \log_2 p(l) - p(r) \log_2 p(r)$$
where $p(l)$ is the proportion of the parent node that belong to the left node, $p(r)$ to the right. The algorithm maximize this entropy, which in our case can be seen as information gain.
Random Forest is an ensemble of decision trees which algorithm (as described in raschka:
If we take the usual linear regression with no constant, we are searching for the weights $W = {w_i^1}, i \in [0,2]$. The variables $X^1 = {X_i^1}, i \in [0,2]$ are called features. We have observations of $X_0^2$ and we want to estimate them from the features $X_i^1, i \in [0,2]$, that is finding $\hat{X}_0^2 = W'X^1$.
To estimate the weights, that is finding $\hat{w}_i^1, i \in [0,2]$, the Ordinary Least Squares method. Our programme is:
With first order conditions, we can find the estimate of $w_i^1, i \in [0,2]$ and we could prove that it is linear, unbiased, and efficient. For simplicity, if there were only one feature, the result would be: $$\hat{w}_0^1= \frac{(X_0^1)' X_0^2}{(X_0^1)'X_0 ^1}$$ which is echoing equation (ref), we could also add a constant in our model.
In his work in perceptron, Frank Rosenblatt proposed in 1957 a perceptron algorithm where for a binary classification task, an activation function $\sigma ( )$ is added in series.
Frank Rosenblatt defined a threshold $\theta$ for which:
This linear classification might not be satisfactory for scattered data sets and we might need to either use more complex activation function or add some layers. In which case we move to Artificial Neural Network.
\quad
We\footnote{perceptron_example.py} propose a simple illustration of the perceptron algorithm. We generated two features, $X_1$ and $X_2$ and we define the output $Y$ as:
That is indeed splitting the "observation" as follow:
With our perceptron rule, we want to find the estimates $\hat{\alpha}_1$, $\hat{\alpha}_2$, and $\hat{\theta}$ (pretending we did not know their values).
\quad We follow the steps:
We propose as an illustrative exercise to implement this algorithm with: $\epsilon = 10^{-4}$ and $\eta = 10^{-2}$ and limit our iteration to $10^4$. In raschka you can find this algorithm in a more pythonic way, defining classes which can be found \href{https://github.com/rasbt/python-machine-learning-book-3rd-edition/blob/master/ch02/ch02.ipynb}{here} and illustrated:
\quad
This separates the data in two categories as follow:
We seperate our data set into a train and test one and compare our prediction accuracy with the scikit-learn Perceptron algorithm one.
\stepcounter{question}
With the perceptron algorithm, we updated observation by observation for each iteration. The threshold activation function we used is not differentiable. The ADAptive LInear NEuron algorithm add a cost function:
as this cost function is convex, we can use the gradient descent algorithm introduced section (ref), we derive the cost function with respect to each weight, we find an update:
You should find this algorithm\footnote{code: adaline_example.py} much faster than the simple perceptron.
In a simple example, we could have $\sigma(x) = \frac{1}{1 + e^{-x}}$ as we use in a logistic regression for the activation function. We would then have: $$X_0^2 = \sigma \left( w_{00}^1 X_0^1 + w_{10}^1 X_1^1 + w_{20}^1 X_2^1 \right) $$
The same philosophy goes for the next layer:
$$X_0^3 = Y_k \left( w_{00}^2 X_0^2 + w_{10}^2 X_1^2 + \ldots \right) $$
We call it deep when there are multiple hidden layers, in order to estimate those models, we use deep learning techniques:
A support vector machine (SVM) constructs a hyper-plane or set of hyper-planes in a high or infinite dimensional space, which can be used for classification, regression or other tasks. Intuitively, a good separation is achieved by the hyper-plane that has the largest distance to the nearest training data points of any class (so-called functional margin), since in general the larger the margin the lower the generalization error of the classifier. We suggest to follow the \href{https://scikit-learn.org/stable/modules/svm.html}{scikit-learn section on SVM to go further}.
We illustrate\footnote{code SVM_illustration.py} an implementation of a SVM on the same data set as what we used for the perceptron example and the logistic regression using raschka typology:
We generate non-linearly separable data: $$ \alpha_1\cos(X_1-.1) + \alpha_2 \sin\left[(1+X_1)X_2 \right] X_2 $$ We then use the radial basis function (rbf) kernel that transform the data set before apply the SVM:
We find the following accuracy scores:
Linear Discriminant Analysis (LDA) is a method similar to Principal Component Analysis (PCA). They reduce the dimensions in a dataset. LDA is unsupervised while PCA is supervised. LDA reduce the dimension by finding feature subspace that optimizes class separability.
We\footnote{PCA_example.py} can resort to a principal component analysis to reduce the dimensions of a data set to a minimum set of common factors. In the case where one factor would explain more than 70% of the variance, we would keep only one factor.
For this exercise, we focus on financial data of the French biotech sector over 5 year since September 2012 with a capitalization over EUR 100 million. We use data directly from Euronext and focus on the following companies: ABIVAX, STALLERGENES GREER, GENFIT, ERYTECH PHARMA, VIRBAC, BOIRON, NICOX REGROUPEMENT, VETOQUINOL, INVENTIVA, VALNEVA, ADOCIA, CELLECTIS, ONXEO, GUERBET, DBV TECHNOLOGIES, GENSIGHT BIOLOGICS, PCAS, INNATE PHARMA, POXEL, PHARNEXT, NANOBIOTIX, AB SCIENCE, TRANSGENE, GENKYOTEX.
We standardize our returns (centre-reduit):
where, for each stock $i$, $\bar{R_i}$ is the average and $\sigma_i$ its standard deviation, and assume the their covariance matrix is $\Sigma$.
\quad
The first principal component is $F_{1,t} = w_1' R_t^{cr}$, where $w_1$ is the eigenvector associated with the largest eigenvalue $\lambda_1$ of $\Sigma$. This value of $w_1$ solves the problem $\max_w w' \Sigma w$ subject to the normalization $w' w = 1$. The eigenvalue $\lambda_1$ equals $VAR(F_{1,t}) = w_1' \Sigma w_1$. The $j$th principal component solves the same problem, but under the additional restriction that $w_i' w_j =0 \forall i< j$. The solution is the eigenvector associated with the $j$th largest eigenvalue $\lambda_j$ (which equals $VAR(F_{j,t}) = w_j' \Sigma w_j$).
We extract factor 1 and 2 ($F_{1,t}$ and $F_{2,t}$) with a principal component analysis :
Then for each stock returns $i$, we use the following linear regression:
We assume that the factors and the residuals are following normal distribution\footnote{This is an area of improvement as during stress event the distribution is not normal, cf. Gourieroux - Dubecq}
We chose a number $S$ of simulations we want to run and follow those 4 steps to generate a shock on the portfolio:
Work on the \href{https://www.eba.europa.eu/documents/10180/1383302/2016+EU-wide+stress+test-Adverse+macro-financial+scenario.pdf}{EBA stress test scenario that were used in 2016}. You have the GDP, inflation and unemployment information stored in \verb|unemployment_EBA_2016.csv|, \verb|GDP_EBA_2016.csv|, and \verb|inflation_EBA_2016.csv|
\quad
Find the most important factor behind this scenarios of 2016. Once available, compare them to the 2018 scenarios.
Boosting\footnote{Text from MLBoE} describes a set of methods which transform, or boost, a weak learner into a strong one (model). This means that one iteratively applies models to the data which individually have weak generalisation properties, but the final ensemble of models generalises well. Two popular boosting techniques are gradient boosting and adaptive boosting, or adaBoost. In gradient boosting, one sequentially updates the estimator in the direction of change of the loss function, i.e. a form of error correction. For adaBoost, one gives larger weights to misclassified or badly fitted observations at each consecutive step until one reaches a certain stopping criterion. Stopping criteria can be linked to target test or a flattening out-of-bag performance. AdaBoost requires a potentially large amount of data to converge and is susceptible to noise in the data. Boosting is often applied to tree models.
We can randomly separate our data set into two sets: a Train and a Test data set. The Train data set will have two third of the data, and the Test set the rest. The practice is then to fit a model on the Train data set, then to use prediction on the test data set and see if the predictions are "close" to our empirical observations.
We refer directly to sklearn tutorial on \href{https://scikit-learn.org/stable/modules/grid_search.html#}{Tuning the hyper-parameters of an estimator} and provide a detailed practical example section (ref)
We sort our results by F1 score.
When lending, we can consider that false negatives (credit worthy clients get rejected) are more acceptable than false positives (non-credit worthy clients get a loan from us).
\resizebox{\textwidth}{!}{
}
We optimize\footnote{code: ML_credit_card_NN.py} the estimator for the neural network which requires computing power. We have to chose:
We find:
Let's imagine you spend some time at Dunkerque end 2023. You want to sail and need to know\footnote{dunkerque_tide.py} the tide's height at any time of the day. You have only access to the time of the low and high tide and both heights.
Hypothesize that the tide height follows:
where $a$ is half the tidal range, $b$ is the mean sea level, $\omega$ is angular frequency, $\phi$ fixes phase, and $t$ the time (in hours) since the low tide. With low tide at $t=0$ and period $T$, you should set $\omega = \tfrac{2\pi}{T}$ and $\phi = -\tfrac{\pi}{2}$ so that $h(0) = b-a$ matches the low tide. Between one low tide and the next high tide, the elapsed time is about 6h12. That corresponds to a quarter of a full tidal period (since full period is 12h25). So $\omega = 2\pi / T$ with $T\approx 12.42$ h. You only have the following information:
This give the tide height evolution on Thursday 28th Dec. 2023:
Note that the tide's coefficient $c$ will impact the tide's height. We import the data set for different tide coefficients and try to see if there is a systematic linear relationship between $a$ and $c$ on one side and $b$ and $c$ on the other.
\stepcounter{question}
To apply the skills developed so far, we suggest to play Liar's Dice and play the \href{http://www.perudo.com/perudo-rules.html}{Perudo}\footnote{code: vansteenberghe_perudo.py} version and take the rules from perudo.com, but with some twists to make it simple:
The next player can also announce "calza" which means that he believes that the announcement of the player just before him was the correct guess. In fact, any player can call "calza" at any time even on his own announcement. If the player calling "calza" is correct, then we give him another dice, but if he is wrong then he loses one of his dice.
Suppose there are 6 players in a game, which means there are initially 30 dice in play (let's denote this number by \( n_{\text{dice}} \)). The goal is to determine the most probable number of dice that show the same face value. This can be achieved through a Monte Carlo simulation where numerous random draws (“tirages”) are performed. The outcome of a single dice roll can range from 1 to 6. However, it's crucial to remember that an ace (represented by the value 1) is special and can stand for any value between 2 and 6, making it analogous to a joker.
From the simulations, we confirm that the average count of similar dice is 10. But what about the probability associated with making a “Dudo” call?
A “Dudo” is correctly called when the actual count of a particular dice face is strictly less than a player's announced count. We can accumulate all such counts that validate a “Dudo” call. For this particular game setup, if your predecessor announces a count of 13, calling a “Dudo” gives you an over 80% probability of being correct.
However, it's essential to recognize that prior players might have provided hints about their dice counts (even though they might be bluffing). This additional information can change the odds.
Now before you speak, several players did before you, but be careful, they might be bluffing. Let's imagine that you have three dices of $2$, one ace and so far two players spoke and told the audience $3$ twos and then $12$ twos. If you say $13$ twos, are the odds still at $80\%$ against you?
\stepcounter{question}
In the following we will study the probability laws governing random variables (rv), their distributions. We are interested in their cumulative distribution function $P(X \leq x)$.
If you really want to understand probability, you'd need to dig to measure theory tao2011 to be able to understand the following definition:
In this book, you'll need to get to the following notions:
For an introduction to probability, rv, see hansen2022.
The probability that a discrete rv $X$ equals the value $x$, $P(X=x)$, is called the probability mass function. The expectation of $X$ is the probability-weighted average. For a discrete rv, taking the value $x_i$ such that $P(X=x_i)=p_i$: $E(X)=\sum_{i=1}^\infty p_i x_i$. For a continuous rv, if $P(X \leq x)$ is differentiable, their probability density function $P(X=x)=\frac{d P(X \leq x)}{dx}$. Notation: we write the pdf $f_X(x)$ and the cdf $F_X(x)$.
We follow kaplan1981quantitative to define frequency versus probability. Frequency is an objective, measurable quantity defined as the outcome of repeated (actual or conceptual) experiments; it belongs to the empirical domain of data and statistics. Probability, by contrast, is a numerical measure of uncertainty or degree of belief, representing a state of knowledge rather than a physical property of the world. Uncertainty is assumed to be comparable across propositions, which allows it to be quantified on a calibrated scale. Frequencies provide the reference standard for this calibration but probability, once defined, applies more broadly, including to unique or non-repeatable events where no frequency data exist. Although probability is subjective in that it reflects information rather than reality, it is objective conditional on information: given the same background knowledge, rational agents must assign the same probabilities. Accordingly, statistics is the science of handling data (frequencies), whereas probability is the science of reasoning under incomplete information. Far from being unusable with limited data, probability is precisely the appropriate tool in such settings.
einhorn1985 writes:
einhorn1985 defines $0\leq \theta \leq 1$ the amount of ambiguity and $\beta$ an individual's attitude toward ambiguity. The judged probability of an initial guessed probability $p_A$ is expressed as:
To explain the Ellsberg's paradow, they suggest taking $\theta>0$ and $\beta<1$
For any data set, as an empiricist, we are making assumptions. The strongest and least credible assumptions we make are on the data-generating process (DGP)\footnote{Davidson1978aa consider that "one could characterise 'econometric modelling' as an attempt to match the hypothetical data generation process postulated by economic theory with the main properties of the observed data."} where we assume that there is an underlying mechanism that is generating the observations under study. This is discusses in breiman2001statistical and we discuss further such approaches with machine learning section (ref).
In our example, this means that we assume that the CAC 40 daily returns are generated by a process, we will introduce possible distributions in the following sections. We further initially assume that the return at date $t$ is drawn from the DGP independently with respect to the outcome at date $t-1$, but for time series we will need later to tackle the fact that this hypothesis doesn't hold.
Further, when applying an ordinary least square method to estimate a linear regression, we will assume (and test) that the residuals have mean $0$, were generated by the same underlying distribution, and are independent of each other. For the last two conditions, we say that the residuals are independent and identically distributed (iid). This is often violated: residuals are serially correlated or their variance differ for some observations (heteroskedasticity).
Armatte2017aa review the work of Edmond Malinvaud and his position on "A Priorism versus Empiricism":
he stresses the
And he gives this recommendation:
We follow section (ref) and use the history of CAC 40 index including dividends (net returns)\footnote{vansteenberghe_samples.py}. Now if as an investor I believe that the history or daily returns are representative of the likely upcoming daily returns, I might want to know to likelihood of a future daily return of $0.045$ (i.e. $4.5\%$). I can approach this question two ways:
In formal terms, I consider the daily returns to be a continuous rv $X$ and believe that I can find a pdf and cdf that is stable over time and can help me predict the realizations.
If we build an empirical histogram, we have to chose the granularity of the bins. If doing it manually, we can start by limiting our representation to 2 bins: $]-\infty,0]$ and $]0,+\infty[$. In our representation, we have:
We can normalize this so that the sum of the bins is equal to $1$ ($100\%$):
In other words, empirically, I have $52\%$ of probability to have a positive daily return\footnote{before concluding, we would need to perform statistical test to see if this is indeed statistically different than $50\%$}. $P(X > 0)=0.52$.
\quad
We can use directly the python command and chose both the number of bins ($20$ in this example) and whether we normalize the histogram or not:
So to answer our initial question on the likelihood to have a daily return of $4.5\%$, we can zoom in:
For a continuous rv $X$, the pdf $P(X=0.045)=0$, hansen2022, we are only able to quantify the probability for $X$ to fall in an interval. We visually observe a probability around $0.3\%$ (so on the y-axis between $0.2$ and $0.4$) which we confirm by counting how many daily return fell inside the bucket $[0.04,0.05]$ of size $0.01$ over the total count of daily return: $0.38\%$. You can look at equation (1.6) in li2023nonparametric.
We will present various parametric distributions in the coming sections\footnote{all of those are unimodal (that is, have a single local maximum), meaning that if you work with a multimodal distributed sample you will need to use other technics e.g. KDE we cover in section (ref)}, we need to believe that it is possible to fit a "function" $P(X=x)$ that gives us the probability that the daily return is $x$. Note that the goodness of fit of that function will be relevant if we want to trust the answer on the probability, for this there exist various tests, one of which is presented section (ref).
We define the expectation as a weighed average: $$E(X) = \int_{-\infty}^\infty x P(X=x) dx$$ or in unified notation $E(X) = \int_{-\infty}^\infty x dF_X(x)$
in our case $E\left(X\right)=0.038\%$.
We define the standard deviation as: $$s(X) = \sqrt{\int_{-\infty}^\infty \left[x-E(X)\right]^2 P(X=x) dx}$$ in our case we find $0.013$
We introduce the normal distribution $\mathcal{N}(\mu,\,\sigma^2)$. Its probability density function, $P(X=x)$, is $$f(x;\mu,\sigma) = \frac{1}{\sqrt{2\pi \sigma^2}} e^{-\frac{(x-\mu)^2}{2\sigma^2}}$$ in simple words this function gives us the probability that a random variable $X$ takes the value $x$ if the random variable follow the normal distribution, we typically write $X \sim \mathcal{N}(\mu,\,\sigma^2)$. For illustration, we take the parameter $\mu$ and $\sigma$ as the empirical mean and standard deviation observed on the CAC 40 daily returns:
\stepcounter{question}
The mean or first moment of the distribution can be computed following Section (ref). We start from
with $f(x;\mu,\sigma)$ the $\mathcal{N}(\mu,\sigma^2)$ density.
Make the change of variable \[ z = \frac{x-\mu}{\sigma}, \qquad x = \mu + \sigma z, \qquad dx = \sigma dz. \]
The integral becomes \[ \mu^1 = \int_{-\infty}^{+\infty} (\mu + \sigma z)\, \frac{1}{\sqrt{2\pi}} e^{-z^2/2}\, dz. \]
Split into two terms: \[ \mu^1 = \mu \underbrace{\int_{-\infty}^{+\infty} \frac{1}{\sqrt{2\pi}} e^{-z^2/2}\, dz}_{=1} + \sigma \underbrace{\int_{-\infty}^{+\infty} z \frac{1}{\sqrt{2\pi}} e^{-z^2/2}\, dz}_{=0}. \]
Hence \[ \mu^1 = \mu. \]
We can draw randomly from a normal distribution (let's say $10^4$ draws) and then plot the histogram against the probability density:
The empirical data doesn't fit the theoretical normal distribution, visually you can observe fatter tails:
Let's imagine that you have a sample of size $n$, $[x_1, \ldots, x_n ]$, taken from a population governed by a normal law. We want to estimate two parameters of the parametric law, $\mu$ and $\sigma$. This means you are searching for the parameters that are making a normal law most likely to have generated the sample you observe. We call this method the Maximum Likelihood Estimation. Said differently:
The probability density function for a sample of $n$ independent identically distributed normal random variables is the likelihood we want to maximize:
as the random variables are independent and identically distributed, equation (ref) can be written:
To say it differently, the probability that the random variable $X$ took the value $x_i$, for any $i$ ($\forall i$) is by definition the probability density function and here we use a parametric density function with two parameters that have to be estimated:
We want to choose $\alpha$ and $\mu$ that maximize $L(\mu,\sigma)$:
As the logarithm function is continuous strictly increasing and it is easier to deal with the log-likelihood, we maximize: $$\log\left( L(\mu,\sigma) \right) = -\frac{n}{2} \log\left( 2 \pi \sigma^2 \right) - \frac{1}{2 \sigma^2}\sum_{i=1}^n(x_i-\mu)^2$$
In a maximization problem, we can apply first order conditions, they yields the following parameters that will maximize our probabilities: $$\hat{\mu} = \sum_{i=1}^n \frac{x_i}{n}$$ $$\hat{\sigma}^2 = \frac{1}{n}\sum_{i=1}^n (x_i-\hat{\mu})^2$$
It can be demonstrated that the estimator of the variance is biased:
Intuition. Why is there a bias? The key point is that in the sample variance formula we subtract the sample mean $\hat{\mu}$, rather than the true mean $\mu$. Because $\hat{\mu}$ itself is estimated from the data, it tends to be closer to the sample points than the true $\mu$. This makes the deviations $(x_i-\hat{\mu})^2$ systematically smaller than the deviations $(x_i-\mu)^2$, and therefore the variance estimator underestimates the true variance on average.
If we divide by $n$ as in the maximum likelihood formula, we get this downward bias of factor $\tfrac{n-1}{n}$. That is why the unbiased sample variance uses $n-1$ in the denominator:
so that \[ E[s^2] = \sigma^2. \]
Interpretation. The difference between dividing by $n$ and $n-1$ becomes negligible for large $n$ (since $\tfrac{n-1}{n} \to 1$ as $n\to\infty$). But for small samples, the correction is important: the unbiased variance estimator avoids the systematic underestimation that comes from using $\hat{\mu}$ in place of $\mu$.
We follow ruppert2015.
The matrix of second partial derivatives of a function is called its Hessian matrix, so the Fisher information matrix is the expectation of the Hessian of the negative log-likelihood:
The standard errors are the square roots of the diagonal entries of the inverse of the Fisher information matrix:
And finally, to be able to estimate a confidence interval:
For a normal law, we have: $$\frac{\partial^2 L}{\partial \mu^2} = -\frac{\sum_i 1}{\sigma^2} = -\frac{n}{\sigma^2}$$ hence $F(\hat{\mu}) = \frac{n}{\sigma^2}$ and $s_{\hat{\mu}} = \frac{\sigma}{\sqrt{n}}$. Thus, $\hat{\mu} \sim \mathcal{N}\left(\mu, \frac{\sigma^2}{n} \right)$ or:
and $$ z_{.05/2} \leq \frac{\hat{\mu} - \mu}{\frac{\sigma}{\sqrt{n}}} \leq z_{1-.05/2} $$ finally
and remember that $z_{.05/2} = - z_{1-.05/2}$.
With our DGP in mind (a normal law for the returns) and our empirical observations, we can answer the previous question for the probability to see a daily return of $4.5\%$, or more specifically, the probability that our CAC 40 daily return will fall between $4$ and $5\%$:
\quad
With a normal distribution we tend to underestimate the probability of having absolute large daily returns.
The cumulative distribution function (cdf) answers the question: "What is the probability that my returns are below or equal to x?"
$$P(X\leq x)=\int_{-\infty}^x P(X) dX$$
We can apply a normal law cdf to the CAC 40 daily return:
With the empirical distribution function, for the cumulative up to $4.5\%$ daily return, we obtain $99.6\%$, but for the normal distribution we find $99.95\%$ where again with the normal distribution we underestimate the probability of having daily returns exceeding $4.5\%$: in practice we had $0.4\%$ of the daily returns above $4.5\%$ but the normal distribution function estimates that probability at $0.05\%$.
\stepcounter{question}
Now let's imagine that we have a random variable $X$ that follows $\mathcal{N}(0,\,\sigma^2)$, we might wonder the mean of the $\vert X \vert$ random variable, $E(\vert X \vert)$. Writing $x=\sigma z$ where conveniently $Z \sim \mathcal{N}(0,\,1)$, we have $E(\vert X \vert) = \sigma E(\vert Z \vert)$. And: $$ E(Z) = \int_{-\infty}^{+\infty} \frac{x}{\sqrt{2\pi}} e^{-\frac{x^2}{2}} dx$$ $$ E(\vert Z \vert) = \frac{2}{\sqrt{2\pi}} \int_{0}^{+\infty} x e^{-\frac{x^2}{2}} dx$$ $$ E(\vert Z \vert) = \sqrt{\frac{2}{\pi}} \left[-e^{-\frac{t^2}{2}} \right]_O^\infty = \sqrt{\frac{2}{\pi}}$$ Finally, $E(\vert X \vert) = \sigma \sqrt{\frac{2}{\pi}}$
We can verify this by drawing a sample and computing the formula to compare.
Let's imagine that we have a random variable $X$ that follows $\mathcal{N}(\mu,\,\sigma^2)$. If you have a sample of size $N$, what is the largest value you expect to observe in this sample?
\stepcounter{question}
Exponential distributions are a particular case of gamma distribution. We write the probability density function of the exponential distribution:
For illustration here the exponential probability density function:
The likelihood function on a sample of size $n$ is: $$ L(x_1,\ldots,x_n\vert \alpha) = \alpha^n \exp\left(-\alpha \sum_{i=1}^n x_i \right) $$
Deriving this expression with respect to $\alpha$ to find an optimum $\alpha_o$: $$n\alpha_o^{n-1} \exp\left(-\alpha_o \sum_{i=1}^n x_i \right) - \sum_{i=1}^n x_i \alpha_o^n \exp\left(-\alpha_o \sum_{i=1}^n x_i \right)=0$$ an optimum is thus find for the reciprocal of the sample mean: $$ \alpha_o = \frac{n}{\sum_{i=1}^n x_i}$$
\quad
For the exponential law the first moment (the mean) is:
and we can demonstrate that the variance of the exponential law is $\frac{1}{\alpha^2}$.
\quad
We can fit this to our negative CAC 40 return (we need to take the absolute values to be able to fit our exponential distribution):
Typically, in insurance claims modelling (non-life), we model the occurrence of claims as a Poisson\footnote{the probability density function of a Poisson distribution is $P(N_t=k) = e^{-\lambda t} \frac{(\lambda t)^k}{k!}, \quad \lambda > 0$} process of parameter $\lambda$ and the magnitudes of claims as an exponential law of parameter $\alpha$. Writing that the claim epoch follow a Poisson process is equivalent to writing that the inter-occurrence time follow an exponential distribution with rate parameter $\lambda$ and a probability density function of $\lambda e^{-\lambda x}$ for $x \geq 0$.
We model an insurance company facing claims arriving over time with a Poisson law $N_t$ and which claim size $C_i$ follow and exponential law. We have the compound process: $X_t=\sum_{i=1}^{N_t} C_i$ and its estimate $\hat{X}_t=\sum_{i=1}^{N_t} \hat{C}_i$.
We have: $E[X(t)] = E[N(t)] E[C] = \frac{\lambda}{\alpha}$ and $Var[X(t)] = E[N(t)] E[C^2] = \lambda ( \sigma^2 + \mu^2 ) = \lambda (\frac{1}{\alpha^4}+\frac{1}{\alpha^2}) $
We introduce the Generalized Extreme Value distribution which has a cumulative distribution function:
for $1+\xi\left(\frac{x-\mu}{\sigma}\right)>0$, $\xi \neq 0$ the shape parameter, $\mu$ the location parameter and $\sigma>0$ the scale parameter. We can derive a probability density function
This distribution is sometimes called the Fisher-Tippett distribution.
For $\xi = 0$, we are back to the Gumbel distribution introduced equation (ref):
We write the probability density function of the Gumbel distribution:
For illustration here the standardized GEV probability density function for different values of $\xi$:
For $\xi \geq 1$ this distribution has neither finite mean nor variance, illustration wih $\xi=1.5$:
To build confidence interval of the GEV estimators, you can follow section (ref) and the expressions of the Fischer information matrix are derived in Prescott:1983ur.
We write the Gamma distribution with distribution function $f(x) = \frac{r^s}{\Gamma(s)}x^{s-1}e^{-r x}$, $x>0$, $s>0$ the shape parameter, $r>0$ the rate parameter and $\Gamma(s) = \int_0^{+\infty}u^{s-1}e^{-u}du$, we have the property\footnote{some elements of demonstration \href{https://en.wikipedia.org/wiki/List_of_integrals_of_exponential_functions}{here}} $\Gamma(s+1)=s \Gamma(s)$ and $\frac{\Gamma(s+1)}{\beta^{s+1}} = \int_0^{+\infty}u^{s}e^{-\beta u}du$. Such distribution can be useful in insurance modelling (claims size distribution). If $s>1$ then it is alike a log-normal distribution:
For a random variable following a gamma distribution, we have: $$E(X) = \int_0^{+\infty} \frac{r^s}{\Gamma(s)} u^{s}e^{-r u}du = \frac{r^s}{\Gamma(s)} \int_0^{+\infty} u^{s} e^{-r u}du$$ $$E(X) = \frac{r^s}{\Gamma(s)} \frac{\Gamma(s+1)}{r^{s+1}} = \frac{s}{r} $$
For the variance, we have $Var(X) = E(X^2) - \left(E(X) \right)^2$ and $E(X^2) = \frac{s (s+1)}{r^2}$ so $Var(X) = \frac{s}{r^2}$.
Let $X$ be a continuous rv with a continuous and invertible (strictly increasing) cdf. Then the random variable $Y=F_X(x)$ has the support $[0,1]$ and $F_Y(y)=y$ and $f_Y(y)=1$, hansen2022.
This means that from a uniform rv, we can generate rv following other distribution. Let's imagine we want to generate a rv following the Pareto cdf:
with $x\geq x_m >0$ and $\alpha>0$. We take\footnote{vansteenberghe_probability_integral_transformation.py} $Y \sim U$ with $U$ the uniform cdf over $[0,1]$. We want $y=1-\left(\frac{x_m}{x} \right)^\alpha$ so we take $x= x_m (1-y)^{-\frac{1}{\alpha}}$.
We start by generating a sample uniformly distributed, so we can check that as expected our expectation is $0.5$ and the standard deviation $.25$, remembering that $P(X=x)=1$ on $[0,1]$ $$E(X) = \int_0^1 x dx = \frac{1}{2}$$ $$s(X) = \sqrt{\int_{0}^1 \left[x- \frac{1}{2} \right]^2 dx} = \sqrt{\left[ \frac{\left( x-\frac{1}{2}\right)^3}{3} \right]_0^1} = \sqrt{\frac{1}{12}}$$
If we have a sample of size $n$ from a population, we consider\footnote{vansteenberghe_Central_Limit_Theorem.py} the random variable $X$ which have the outcome values $x$ (that would be $x \in [0,\infty [$ for a stock price), each outcome having the probability of occurring of $p(x)$. We can consider that this random variables are taken from a greater population with mean $\mu$ and standard deviation $\sigma$ (although a price process would never follow a normal law). In practice, we do not have access to the full population and have to work with finite samples.
Not that the outcome values $X$ are considered a random sample of independent and identically distributed random variables with the same probability distribution function. We will see in section (ref) possible distributions function for financial returns.
In practical terms, if over 5 days we observe the following price values $x_t$ for a stock $X$, we have the probabilities:
This is a very basic example and in practice you would need much more observations to work with plausible probability values.
The expected value of your random variable $X$ is called the sample mean, we can show that it is unbiased compared with the population true mean:
We can also write:
In our example, this would be: $E(X) = 1.1 \times \frac{1}{5} + 1.2 \times \frac{3}{5} + 1.3 \times \frac{1}{5} = \frac{6}{5} = 1.2$
We can state the weak law of large numbers:
Statement: Let $X_1, X_2, \ldots$ be a sequence of i.i.d. random variables with expected value $E(X_i) = \mu$ and variance $Var(X_i) = \sigma^2 < \infty$. Then, for any $\epsilon > 0$, \[ \lim_{n \to \infty} P\left( \left| \frac{1}{n} \sum_{i=1}^{n} X_i - \mu \right| > \epsilon \right) = 0 \] Illustration: Imagine flipping a fair coin, where Heads is `1` and Tails is `0`. The expected value is $ \mu = 0.5 $. If you flip the coin many times and take the average of the results, this average will converge to 0.5 as the number of flips goes to infinity.
Another more mathematical way of writing the WLLN: If we have $X_i, i\in[1,n]$ iid, with $E\left( \vert X \vert \right)<\infty$, then $\bar{X}_n = \frac{1}{n} \sum_{i=1}^n X_i \overset{p}{\to} E(X)$.
\quad
If we further assume $E(X)=\mu$ and $\sigma(X) = \sigma$. The variance of the mean of our sample is biased:
If we want an unbiased standard deviation\footnote{this can be measured using the standard deviation function with a degree of freedom of 1 ($\mbox{ddof} = 1$)} measure from our sample, one can demonstrate that it can be written as:
The variance of the sample is:
In our example, this would be: $V(X) = (1.1-1.2)^2 \times \frac{1}{5} + (1.2-1.2)^2 \times \frac{3}{5} + (1.3-1.2)^2 \times \frac{1}{5} = 0.004$
Statement: Let $X_1, X_2, \ldots$ be i.i.d. random variables with $E(X_i) = \mu$ and $Var(X_i) = \sigma^2 < \infty$. Then, as $n \to \infty$, the random variable \[ Z_n = \frac{\sum_{i=1}^{n} X_i - n\mu}{\sigma \sqrt{n}} \] converges in distribution to the standard normal distribution, i.e., \[ Z_n \xrightarrow{d} Z \sim N(0,1) \]
\stepcounter{question}
A more compact way to write the CLT: If we have $X_i, i\in[1,n]$ iid, with $E\left( \vert X \vert \right)<\infty$, with $E(X) = \mu$ and $\sigma(X) = \sigma$, then if we define $\bar{Z}_n = \frac{\bar{X}_n - \mu}{\sqrt{\frac{\sigma^2}{n}}}$, then $\bar{Z}_n$ converges in distribution, $n \to \infty$ to the standard normal distribution $\mathcal{N}(0,\,1)$.
\quad
We suggest to illustrate this with samples taken from a population that has a uniform distribution function, with its probability distribution function:
we take $a=0$ and $b=1$.
\stepcounter{question}
Gosset has established under the pseudonym Student that:
$$\frac{\bar{X}_n-\mu}{\sqrt{\frac{\sigma^2}{n}}} \sim \mathcal{N}(0,\,1)$$
and $t_{n-1}$ being the t-Student distribution with $n-1$ degrees of freedom:
$$\frac{\bar{X}_n-\mu}{\sqrt{\frac{S_n^2}{n}}} \sim t_{n-1}$$
And remember that for a t-Student distribution with $\nu$ degrees of freedom $t_\nu$, then when $\nu \to \infty$, then it becomes a standard normal distribution $\mathcal{N}(0,\,1)$:
This means that the greater our sample size $n$, the closer to the actual population mean and standard deviation we will measure. This can be rewritten and is called the Central Limit Theroem:
The standard deviation of our sample is defined as:
its skewness is defined as:
its kurtosis is defined as:
The covariance between two random variables $X$ and $Y$ indicates how likely they are to occur together:
If we look at an observation vector $X=[x_0,x_1,\ldots,x_{T-1}]$ (like the close price of a company's stock), we can call $E(X)$ the expected price, which in practice will be the average over the sample: $$E(X)=\frac{1}{T}\sum_{i=0}^{T-1} x_i$$
The standard deviation is a measure of risk expressed in $\%$: $$\sigma(X)=\sqrt{\frac{1}{T}\sum_{i=0}^{T-1} \left(x_i-E(X) \right)^2}$$
The covariance of $X$ with another vector of observations $Y=[y_0,y_1,\ldots,y_{T-1}]$ can be computed as: $$cov(X,Y)=\frac{1}{T}\sum_{i=0}^{T-1} \left[x_i-E(X) \right] \left[y_i-E(Y) \right]$$
The Pearson coefficient of correlation between $X$ and $Y$ is then: $$corr(X,Y)=\frac{cov(X,Y)}{\sigma(X) \sigma(Y)}$$
We can detail the correlation measures:
with:
One might wonder: should correlation be applied to prices or returns? You can find detailed advice to apply correlation on asset returns \href{https://quantdare.com/correlation-prices-returns/}{here}, \href{http://t.www.na-businesspress.com/JBD/HaberJ_Web12_3_.pdf}{here} and \href{https://www.portfolioprobe.com/2011/01/12/the-number-1-novice-quant-mistake/}{here}.
I follow here johnson1987multivariate.
We have seen how to generate random vectors from a uniform distribution. Now we want to generate a variable $X$ whose cumulative distribution function is specified by $F$.
The inverse distribution function method has two steps:
In the event that $F$ is strictly increasing, then $F^{-1}$ is well defined and step 2 becomes: $X=F^{-1}(U)$.
I test\footnote{numpy_exercise_part4.py} the following two:
When no closed form of $F^{-1}$ is available, the following recursion\footnote{Cf the Newton's method equation (ref).} frequently converges to a fixed point $x^*$: $$x_{n+1}=x_n - \frac{F(x_n)-U}{f(x_n)}$$
A standard normal law can be generated with two independent uniform vectors $U_1$ and $U_2$: $$X = \sqrt{-2 \ln(U_1)}\cos(2\pi U_2)$$ or $$X = \sqrt{-2 \ln(U_1)}\sin(2\pi U_2)$$
Once a standard normal variable $X$ is generated, a normal variable $Y\sim \mathcal{N}(\mu,\sigma^2)$ can be generated with $Y = \sigma X + \mu$
Finally, a log-normal variable can be generated from a normal variable taking its exponential.
A $p$-dimensional random vector $X=\left(X_1,\dots,X_p \right)'$ is defined to have a multivariate normal distribtuion if and only if every nontrivial linear combination of the $p$-components of $X$ has a univariate normal distribution. The distribution of $X$ is denoted $\mathcal{N}_p\left(\mu, \Sigma \right)$, with $\mu_i=E(X_i)$ and $\Sigma_{i,j}=cov(X_i,X_j)$. If $\Sigma$ has full rank $p$, then the density function of $X$ exists. Let $L$ be the Choleski factorization which is then unique lower triangular such that $LL'=\Sigma$. Take $p$ univariate independent random vector $Z=\left(Z_1,\dots,Z_p \right)'$ with $Z\sim \mathcal{N}_p(0,I)$ if I define $X=LZ+\mu$, then $X\sim \mathcal{N}_p\left(\mu, \Sigma \right)$.
As shown in cario1997modeling and henderson2000generating, it is not easy to determine the required $\Sigma_Z$ to get $\Sigma_X$.
The NORTA generation procedure works as follow:
with $F_{X_i}$ the desired cdf, and indeed, the central problem is to select $\Sigma_Z$ that gives $\Sigma_X$ and for most applications there is no closed form: $$\Sigma_X^{ij}=cov(X_i,X_j)=cov\left(F^{-1}_{X_i}\left[\Phi\left( Z_i \right) \right],F^{-1}_{X_j}\left[\Phi\left( Z_j \right) \right] \right)$$
We aim to recover the subjective probability distribution of average inflation over $n$ years. Let \[ \pi(n)\quad\text{be the $n$-year average annual CPI inflation,}\qquad f^{\mathbb{Q}^{n}}_{\pi(n)}(\cdot)\quad\text{its pdf under the $n$-forward measure}. \] There are two complementary routes. First, use zero-coupon inflation caps/floors across strikes to extract the risk-neutral (forward-measure) distribution by differentiating option prices with respect to strike (Breeden--Litzenberger) and discounting with the $n$-year zero rate breeden1978prices,ross1976options,kitsul2013economics. Second, use survey data to obtain a physical or stated-belief distribution and compare it to the options-implied one.
To extract a probability distribution of future inflation from option prices, remember three key distinctions:
\quad
kitsul2013economics define a zero-coupon inflation cap as a contract entered into at time $t$. The seller of the cap promises to pay a fraction \[ \max\!\left( (1+\pi(n))^n - (1+k)^n,\, 0 \right) \] of a notional underlying principal as a single payment in $n$ years' time, where $\pi(n)$ denotes the average annual total consumer price index (CPI) inflation rate from $t$ to $t+n$ (not seasonally adjusted) and $k$ denotes the strike of the cap. Without loss of generality, the notional underlying principal is normalized to $1$. In exchange for this, the buyer makes an up-front payment of $P_t(k,n)$. The contract is effectively a call option on inflation. If $\pi(n)<k$, then the option is out-of-the-money.
Consider a one-year zero-coupon inflation cap with strike $k=2\%$. If realized inflation is $\pi(1)=3\%$, the payoff is \[ (1+\pi(1))-(1+k)=1.03-1.02=0.01, \] that is, one cent per dollar of notional. If both parties are certain that $\pi(1)=3\%$ (i.e.\ $\mathbb{P}[\pi(1)=3\%]=1$) and there is no risk premium, then the fair price equals the discounted expected payoff: \[ P_t(0.02,1)=e^{-y(1)}\times 0.01. \]
Principle. A price equals the expected present value of its future payoff under a pricing measure. If $y(n)$ is the continuously compounded $n$-year zero rate observed at $t=0$ (so here we take a bond with known return at time $t$ up to maturity and not the risk-free rate, which is unknown at time $t$ up to $t+n$), the price of one unit delivered at $t=n$ is \[ B(0,n)=e^{-n y(n)}. \] This follows from continuous compounding: investing $e^{-n y(n)}$ today ($t=0$) at rate $y(n)$ for $n$ years grows deterministically to $1$ at $t=n$.
Suppose I receive $m$ interest payments at a rate of $\frac{r}{m}$ per annum. After $t$ years, I will receive: $$ \left( 1 + \frac{r}{m} \right)^{mt}$$ which can be rewritten as $\exp\left(m t \log\left[ 1 + \frac{r}{m} \right] \right) \simeq \exp\left( rt \right)$, when $ \frac{r}{m}$ is in the vicinity of 0.
We can get the same approximation with a Taylor expansion. If you have the amount $M(t)$ and invest it at interest rate $r$, then $$M(t+dt) = M(t) + \frac{d M(t)}{dt} \frac{t + dt - t}{1!} + \ldots$$ hence $M(t+dt) - M(t) = M(t) r dt \simeq \frac{d M(t)}{dt} dt$, and $$\frac{d M(t)}{dt} = r M(t)$$ which has the solution $M(0)e^{rt}$.
Applied to the cap. The zero-coupon inflation cap pays only at $t+n$. Therefore
For $n=1$, the payoff is $(\pi(1)-k)_{+}$, so
Integral form. Let $f^{\mathbb{Q}^{1}}_{\pi(1)}$ be the pdf of $\pi(1)$ under $\mathbb{Q}^{1}$. Then
\paragraph{Meaning of $\mathbb{Q}^{1}\!\big(\pi(1)>k\big)$.} $\mathbb{Q}^{1}$ denotes the $1$-forward (or risk-neutral) probability measure. Thus $\mathbb{Q}^{1}\!\big(\pi(1)>k\big)=1-F^{\mathbb{Q}^{1}}_{\pi(1)}(k)$ is the tail probability of $\pi(1)$ under $\mathbb{Q}^{1}$.
\paragraph{Differentiation with a $k$-dependent lower limit.} From the integral form \[ P_t(k,1)=e^{-y(1)}\int_{k}^{\infty}(\pi-k)\,f^{\mathbb{Q}^{1}}_{\pi(1)}(\pi)\,\mathrm{d}\pi, \] set \[ g(k):=\int_{k}^{\infty}(\pi-k)\,f^{\mathbb{Q}^{1}}_{\pi(1)}(\pi)\,\mathrm{d}\pi. \] Apply Leibniz’s rule with $h(\pi,k):=(\pi-k)f^{\mathbb{Q}^{1}}_{\pi(1)}(\pi)$, $a(k)=k$, and $b(k)=\infty$: \[ \frac{\mathrm{d}}{\mathrm{d}k}\int_{a(k)}^{b(k)} h(\pi,k)\,\mathrm{d}\pi = h\big(b(k),k\big)b'(k)-h\big(a(k),k\big)a'(k)+\int_{a(k)}^{b(k)}\frac{\partial h}{\partial k}(\pi,k)\,\mathrm{d}\pi. \] Here $b'(k)=0$, $a'(k)=1$, and $h\big(a(k),k\big)=(k-k)\,f^{\mathbb{Q}^{1}}_{\pi(1)}(k)=0$, while \[ \frac{\partial h}{\partial k}(\pi,k)=-\,f^{\mathbb{Q}^{1}}_{\pi(1)}(\pi). \]
Here, we have \[ \frac{\partial P_t(k,1)}{\partial k} =-\,e^{-y(1)}\,\mathbb{Q}^{1}\!\big(\pi(1)>k\big), \qquad \frac{\partial^{2} P_t(k,1)}{\partial k^{2}} =e^{-y(1)}\,f^{\mathbb{Q}^{1}}_{\pi(1)}(k). \]
\[ \qquad\Rightarrow\qquad \boxed{ f_{\pi(1)}^{\mathbb{Q}^1}(k) = e^{y(1)} \frac{\partial^2 P_t(k,1)}{\partial k^2} }. \]
kitsul2013economics write:
\paragraph{Step 0: Notation.} Let $\pi(n)$ be the average annual CPI inflation from $t$ to $t+n$. For a zero-coupon $n$-year inflation cap with strike $k$, define \[ X \equiv (1+\pi(n))^{n}, \qquad K \equiv (1+k)^{n}. \] The cap payoff at $t+n$ is $(X-K)_{+}$.
\paragraph{Step 1: Pricing under the forward measure (foundation in ross1976options).} Let $y(n)$ be the continuously compounded $n$-year zero yield and $\mathbb{Q}^{n}$ the $n$-forward measure. The time-t price equals discounted forward-measure expectation:
\paragraph{Step 2: Breeden--Litzenberger identity on strikes breeden1978prices,ross1976options.} View $C(K)\equiv e^{-n y(n)}\,\mathbb{E}^{\mathbb{Q}^{n}}[(X-K)_{+}]$ as a call on $X$ with strike $K$. Then
where $f^{\mathbb{Q}^{n}}_{X}$ is the pdf of $X$ under $\mathbb{Q}^{n}$.
\paragraph{Step 3: Chain rule to work in inflation strike $k$.} Since $K=(1+k)^{n}$, define $P(k)\equiv P_t(k,n)=C(K(k))$. By the chain rule,
\paragraph{Step 4: Nonparametric recovery method li2009nonparametric.} Given option quotes $\{(k_i,y_i)\}_{i=1}^{L}$ at one maturity and date, recover a smooth function $k \mapsto P_t(k,n)$ whose second derivative yields the pdf via (ref). The Li--Zhao framework estimates $P$ under slope/convexity constraints implied by (ref)--(ref). Higher convexity of $P(k)$ at a given strike means markets assign more probability mass near that strike. Left-skewed pdf $\Rightarrow$ greater deflation concern.
\paragraph{Step 5: Local polynomial regression to estimate derivatives ait2003nonparametric.} Estimate $P(k)$ and its derivatives by local polynomial regression. Around a target $k$, fit a quadratic locally:
with kernel $K$ and bandwidth $h$. Then \[ \widehat{P}(k)=\hat\beta_0(k), \qquad \widehat{P'}(k)=\hat\beta_1(k), \qquad \widehat{P''}(k)=\hat\beta'_1(k), \]
\paragraph{Step 6: Putting it together with discounting gurkaynak2007us.} Obtain $y(n)$ from the Treasury zero-coupon term structure. Then, for $n=1$,
For $n>1$, combine (ref) with the change of variables \[ f^{\mathbb{Q}^{n}}_{\pi(n)}(\pi) = n(1+\pi)^{n-1}\, f^{\mathbb{Q}^{n}}_{X}\!\big((1+\pi)^{n}\big). \]
\paragraph{No-arbitrage filters used in practice.} Impose $P'(k)\le 0$ (monotonicity in strike) and $P''(k)\ge 0$ (convexity) before derivative extraction. Replace observed quotes by nearby values that satisfy these constraints, then re-estimate via (ref).
\paragraph{Main takeaway for students.} Recover the pdf from option quotes in three moves: (i) price as discounted forward-measure expectation (ref); (ii) differentiate in strike (Breeden--Litzenberger) (ref) and map to inflation strike via the chain rule (ref); (iii) estimate the required derivatives nonparametrically with local polynomials (ref), using $y(n)$ from the term structure.
Let's consider $H$ (height) and $W$ (weight) two continuous rv, the joint distribution function of $(H,W)$ is $F(h,w)=P\left(H \leq h, W \leq w \right)$. The rv $H$ and $W$ are statistically independent if
If we believe that height and weight are statistically independent, then we should\footnote{code: vansteenberghe_height_normal_distrib.py} just study both rv separately. We apply a Chi-square test of independence of variables and reject $H_0$: "there is no significant relationship between weight and height".
In the following, we will try to estimate a relationship between height and weight.
\quad
Hermanussen2001aa found that
We further explore the concept of DGP introduced secion (ref). Now, the weight and height of an individual is likely to be related to his parent's, but appart from those relations, we are closer to an i.i.d. process than with a time series.
In the following, we use Indian 2005-6 Demographic and Health Surveys and would like to thank Professor Dean Spears for his help on collecting the data:
\quad
\quad
Our sample body weights are not normally distributed, but skewed to the right:
We apply a tow-sided t-test and we do not reject the null hypothesis H0: the expected value of this sample made of (presumably) independent observations is equal to $\mu=567.5$.
We apply a two-sided Kolmogorov-Smirnov (cf. section (ref), alternatively we can use the Anderson-Darling test) test on our sample with the null hypothesis: $H_0$ the sample was taken from a normal distributed population and we reject the null hypothesis.
The Lilliefors test is a modification of the Kolmogorov–Smirnov (KS) test, designed to test whether a sample comes from a normal distribution when the mean and variance are unknown and must be estimated from the data.
\paragraph{Step 1. Null hypothesis.} We want to test \[ H_0: X_1, \ldots, X_n \sim \mathcal{N}(\mu,\sigma^2) \] for some unknown parameters $\mu$ and $\sigma^2$.
\paragraph{Step 2. Standardization.} From the data we compute the sample mean and variance: \[ \hat{\mu} = \frac{1}{n}\sum_{i=1}^n X_i, \qquad \hat{\sigma}^2 = \frac{1}{n-1}\sum_{i=1}^n (X_i - \hat{\mu})^2. \] We then standardize: \[ Z_i = \frac{X_i - \hat{\mu}}{\hat{\sigma}}. \]
\paragraph{Step 3. Empirical distribution.} Order the standardized data: \[ Z_{(1)} \leq Z_{(2)} \leq \cdots \leq Z_{(n)}. \] Define the empirical cumulative distribution function (ECDF): \[ F_n(z) = \frac{1}{n} \sum_{i=1}^n \mathbf{1}\{Z_i \leq z\}. \]
\paragraph{Step 4. Compare with fitted normal CDF.} Let $\Phi(z)$ denote the CDF of the standard normal $\mathcal{N}(0,1)$. Compute the Kolmogorov–Smirnov type statistic: \[ D_n = \sup_{z \in \mathbb{R}} \, \big| F_n(z) - \Phi(z) \big|. \]
\paragraph{Step 5. Distribution under $H_0$.} If $\mu$ and $\sigma$ were known, $D_n$ would follow the classical KS distribution. But here they are estimated from the sample, which makes the critical values different. Lilliefors (1967) computed these critical values by Monte Carlo simulation.
\paragraph{Step 6. Decision rule.} Reject $H_0$ at significance level $\alpha$ if \[ D_n > c_\alpha(n), \] where $c_\alpha(n)$ is the Lilliefors critical value.
If taking the log of the weight seems closer to a normal distribution, it still fails the KS test:
Taking the log also reduces the weights skewness.
\quad
Our sample body height is almost normally distributed:
We might want to ask ourselves the question: Are the weights normally distributed? For the normal law, this is so common that several tests have been developped.
We can test normality of the distribution with a Jarque-Bera test, with the test statistics that follows asymptotically a $\chi_2^2$: $$JB = \frac{T}{6} skewness^2 + \frac{T}{24} (kurtosis - 3)^2$$
The statistics is above the critical $\chi_2^2$ threshold and we reject $H_0$ the null hypothesis of normal distribution.
The chi-square test tests if a sample of data came from a population with a specific distribution. The null hypothesis of the test is $H_0$: the data follow a specified distribution.
For the chi-square goodness-of-fit computation, the data are divided into k bins and the test statistic is defined as:
where $O_i$ is the observed frequency for bin $i$ and $E_i$ is the expected frequency for bin $i$.
We reject $H_0$, and visually:
Pearson's chi-squared test is used to test for independence between categorical variables. It compares the observed frequencies of events to the frequencies expected under the null hypothesis. In the context of testing whether a coin is fair, we can use this test to check whether the observed frequencies of heads and tails are consistent with the expected frequencies of a fair coin.
The null hypothesis is that the coin is fair, meaning that the probability of heads ($p_H$) and tails ($p_T$) are both $0.5$. This leads to the following hypotheses:
\[ H_0: p_H = p_T = 0.5 \] \[ H_1: p_H \neq p_T \]
To perform the Pearson's chi-squared test, we conduct a series of coin flips and observe the number of heads and tails. Let $O_H$ and $O_T$ represent the observed counts of heads and tails, respectively, and let $E_H$ and $E_T$ be the expected counts under the null hypothesis. The expected counts are calculated based on the assumption that the coin is fair. For a total of $n$ flips, we expect:
\[ E_H = E_T = \frac{n}{2} \]
The test statistic is then calculated as:
\[ \chi^2 = \frac{(O_H - E_H)^2}{E_H} + \frac{(O_T - E_T)^2}{E_T} \]
This statistic follows a chi-squared distribution with 1 degree of freedom. If the test statistic exceeds the critical value from the chi-squared distribution for a chosen significance level (e.g., $\alpha = 0.05$), we reject the null hypothesis and conclude that the coin is not fair.
In practice, we would compute this test after conducting the experiment and then compare the result to a chi-squared critical value table.
While both the chi-square goodness-of-fit test and Pearson's chi-squared test are used to compare observed and expected frequencies, the goodness-of-fit test is used to check whether data follows a specific distribution, whereas Pearson's chi-squared test is used to assess the independence between categorical variables.
Pearson's chi-squared test can also be applied in the context of classification problems, where the objective is to test whether the frequencies of true positives, false positives, true negatives, and false negatives deviate from what is expected under the assumption that the classifier is random or unbiased. This is particularly useful in evaluating the performance of binary classifiers in fields such as medical diagnostics, where true and false positives matter greatly.
In this context, the observed data is presented in a $2 \times 2$ contingency table, which shows the counts for the following outcomes:
The null hypothesis in this case is that the predictions are independent of the actual classification, meaning that the proportions of true positives, false positives, true negatives, and false negatives do not deviate from what would be expected under a random or unbiased classification. This leads to the hypotheses:
\[ H_0: \text{Predicted classification is independent of actual classification} \] \[ H_1: \text{Predicted classification is not independent of actual classification} \]
To perform the Pearson's chi-squared test, we calculate the expected counts for each of the four cells in the contingency table under the assumption of independence. Let $O_{ij}$ represent the observed frequency for cell $(i,j)$ and $E_{ij}$ the expected frequency for the same cell. The expected frequency for each cell is calculated as:
\[ E_{ij} = \frac{(\text{Row Total}_i)(\text{Column Total}_j)}{\text{Total Observations}} \]
The test statistic is then computed as:
\[ \chi^2 = \sum_{i=1}^{2} \sum_{j=1}^{2} \frac{(O_{ij} - E_{ij})^2}{E_{ij}} \]
This statistic follows a chi-squared distribution with 1 degree of freedom. If the test statistic exceeds the critical value from the chi-squared distribution for a given significance level, we reject the null hypothesis and conclude that there is a significant relationship between predicted and actual classifications.
In the context of true and false positives, this test helps determine whether the classifier's performance is better than random chance or whether there is a bias in the predictions.
The Kolmogorov-Smirnov test allows us to compare the distribution of two observed returns series.
We can implement a two-sample Kolmogorov-Smirnov test, where the null hypothesis is that both groups were sampled from populations with identical distributions.
The two-sided test uses the maximum absolute (vertical) difference between the cumulative distribution function (CDF, denoted $\hat{F}_1$ and $\hat{F}_2$ here) of the distributions of the two data vectors. The test statistic is: $$KS=\max_x \lvert \hat{F}_1(x) - \hat{F}_2(x) \lvert $$
If the K-S statistic is small (but be cautious in how you define small as we talk about differences in distribution functions) or the p-value is high, then we cannot reject the hypothesis $H_0$ that the distributions of the two samples are the same.
\quad
We reject the hypothesis that the daily mean return of the CAC 40 index are normally distributed:
Note that when we apply a two-sample test with a randomly generated normal sample, we do not systematically reject the null hypothesis at the $1\%$ threshold that both sample were taken from the same distribution\footnote{we later move on to the Lilliefors test}.
Definition: A Q-Q plot is a graphical tool used to compare the quantiles of two distributions. Typically, one of the distributions is empirical (from the data), and the other is theoretical. By plotting the quantiles of the empirical distribution against the quantiles of the theoretical distribution, we can visually assess how well the data conforms to the theoretical distribution. If the data lies along a straight line (usually a 45-degree line), it suggests that the data is well-represented by the theoretical distribution.
Illustration: Imagine we have a dataset and want to check if it follows a parametric distribution. We can:
If the points in the Q-Q plot lie roughly along a straight line, it suggests that the data generating process could be the parametric distribution.
Quantile-Quantile plots\footnote{vansteenberghe_qqplots.py} are an important tool to select a theoretical distribution to apply on a given sample. It has both the advantage and the drawback to be a visual tool.
In our approach, we assume that we are observing realisations $x_t$ from a random variable $X$ with theoretical probability distribution function $F$, a quantile of order $p$ is defined as\footnote{when the random variable $X$ is continuous, then $Pr(X<x)=Pr(X \leq x)$} $$Q(p) = F^{-1}(p) = \min {x \in \mathbb{R} : F(x) \geq p} = \min {x \in \mathbb{R} : Pr(X<x) \geq p}$$ If we have $n$ observations, we can write:
\quad
\quad
If there are $i,j$ such that $x_i = x_j$, then we regroup them into $y_1,\ldots,y_k$ for $k\leq n$, we now have, with $n_j$ the number of $x$'s regrouped in $y_j$:
\quad
\quad
We define the empirical distribution function (EDF) as \[ \hat{F}_n(x) \;=\; \frac{1}{n}\sum_{i=1}^n \mathbf{1}_{\{X_i \leq x\}}, \] where $\mathbf{1}_{\{X_i \leq x\}}$ denotes the indicator function.
According to the law of large number\footnote{cf. section (ref)}, this empirical function converges to the theoretical repartition function in probability.
We can now define the empirical quantile of order $p$: $$\hat{Q}_n(p) =\min \{x \in \mathbb{R} : \hat{F}_n(x) \geq p\}$$
$\hat{Q}_n(p)$ converges in probability to $Q(p)$.
Given an ordered sample of independent observations:
A quantile-quantile plot consists of the points:
If we take $n$ observations from a standard normal population, we can order those observations (we should also regroup the observations of same values). Then for each ordered observation $i$, we can compute the quantile of the theoretical standard normal distribution of probability $\frac{i}{n}$.
We run\footnote{code: vansteenberghe_height_normal_distrib.py} a loop searching for alternative distributions and based on the p-values of Kolmogorov-Smirnov test, this can be very specific to the time frame we look at and might not be generalized, we select:
We study here\footnote{vansteenberghe_weight_height_OLS.py} linear regression estimated by Ordinary Least Squares (OLS).
Definition. Linear regression models the relationship between a dependent variable and one or more explanatory variables. In the simple linear case, the model writes \[ y_i = \beta_0 + \beta_1 x_i + \varepsilon_i, \qquad i=1,\dots,n, \] where:
In matrix form: \[ \mathbf{y} = \mathbf{X}\boldsymbol{\beta} + \boldsymbol{\varepsilon}. \]
In Python or R, the formula notation is:
which implicitly includes an intercept. To estimate a regression without an intercept:
Omitting the intercept imposes the restriction $E(y|x=0)=0$. If this restriction is false, the model is misspecified and slope estimates may be biased. It is generally unclear why an economist would estimate a regression without including an intercept. In academic papers, the intercept is sometimes omitted from the tables for brevity, but it was nevertheless included in the actual estimation. In such cases no bias arises, since the intercept is estimated but simply not reported because it is not relevant for the specific point being discussed.
We state the standard assumptions underlying OLS in the cross-sectional case.
The model is linear in the coefficients: \[ y_i = \beta_0 + \beta_1 x_i + \varepsilon_i. \] Nonlinear transformations of $x_i$ are allowed (e.g.\ $x_i^2$, $\log x_i$), provided the model remains linear in $\boldsymbol{\beta}$.
$\{(y_i,x_i)\}_{i=1}^n$ is a random sample from the population.
The regressors are not perfectly linearly dependent. In the simple regression case, this requires \[ \sum_{i=1}^n (x_i-\bar{x})^2 > 0. \]
\[ E(\varepsilon_i \mid x_i) = 0. \] This is the key assumption ensuring unbiasedness of OLS. It implies \[ E(\hat{\beta}_1) = \beta_1. \] It is not tested by a $t$-test on residuals; rather, it is a structural assumption about the data-generating process.
\[ \operatorname{Var}(\varepsilon_i \mid x_i) = \sigma^2. \] Under homoskedasticity, OLS is the Best Linear Unbiased Estimator (BLUE). A common diagnostic is the Breusch--Pagan test based on \[ \hat{\varepsilon}_i^2 = \alpha + \gamma x_i + \nu_i, \] with null hypothesis $H_0: \gamma = 0$.
\[ \operatorname{Cov}(\varepsilon_i,\varepsilon_j \mid X) = 0, \quad i \neq j. \] In time series settings, this can be assessed using Durbin--Watson or Breusch--Godfrey tests.
For exact $t$ and $F$ inference in small samples, one may additionally assume \[ \varepsilon_i \sim \mathcal{N}(0,\sigma^2). \] Normality is not required for unbiasedness or consistency, but simplifies hypothesis testing in finite samples.
Under (A1)--(A4):
Under (A1)--(A5), OLS is BLUE (Gauss--Markov theorem). With normality, exact finite-sample inference based on $t$ and $F$ statistics is valid.
Consider the simple linear regression model \[ y_i=\beta_0+\beta_1 x_i+\varepsilon_i,\qquad i=1,\dots,n, \] and define the OLS estimators $(\hat\beta_0,\hat\beta_1)$ as the minimizers of the sum of squared residuals: \[ (\hat\beta_0,\hat\beta_1) =\arg\min_{\beta_0,\beta_1}\; S(\beta_0,\beta_1), \qquad S(\beta_0,\beta_1)=\sum_{i=1}^n \bigl(y_i-\beta_0-\beta_1 x_i\bigr)^2. \]
Differentiate $S(\beta_0,\beta_1)$ with respect to $\beta_0$ and $\beta_1$.
\paragraph{FOC w.r.t.\ $\beta_0$.} \[ \frac{\partial S}{\partial \beta_0} = \sum_{i=1}^n 2\bigl(y_i-\beta_0-\beta_1 x_i\bigr)(-1) = -2\sum_{i=1}^n \bigl(y_i-\beta_0-\beta_1 x_i\bigr)=0. \] Equivalently, \[ \sum_{i=1}^n y_i - n\beta_0 - \beta_1\sum_{i=1}^n x_i=0 \quad\Longrightarrow\quad \hat\beta_0=\bar y-\hat\beta_1\bar x, \] where $\bar y=\frac{1}{n}\sum_{i=1}^n y_i$ and $\bar x=\frac{1}{n}\sum_{i=1}^n x_i$.
\paragraph{FOC w.r.t.\ $\beta_1$.} \[ \frac{\partial S}{\partial \beta_1} = \sum_{i=1}^n 2\bigl(y_i-\beta_0-\beta_1 x_i\bigr)(-x_i) = -2\sum_{i=1}^n x_i\bigl(y_i-\beta_0-\beta_1 x_i\bigr)=0. \]
Substitute $\hat\beta_0=\bar y-\hat\beta_1\bar x$ into the second FOC: \[ \sum_{i=1}^n x_i\Bigl(y_i-(\bar y-\hat\beta_1\bar x)-\hat\beta_1 x_i\Bigr)=0. \] Use $\sum_{i=1}^n x_i(y_i-\bar y)=\sum_{i=1}^n (x_i-\bar x)(y_i-\bar y)$ (because $\sum_{i=1}^n (y_i-\bar y)=0$), and note that \[ \sum_{i=1}^n x_i(x_i-\bar x)=\sum_{i=1}^n (x_i-\bar x)^2. \] Therefore, \[ \sum_{i=1}^n (x_i-\bar x)(y_i-\bar y) - \hat\beta_1\sum_{i=1}^n (x_i-\bar x)^2=0, \] so that \[ \hat\beta_1 = \frac{\sum_{i=1}^n (x_i-\bar x)(y_i-\bar y)} {\sum_{i=1}^n (x_i-\bar x)^2}. \]
Start from the model and its sample mean: \[ y_i=\beta_0+\beta_1 x_i+\varepsilon_i, \qquad \bar y=\beta_0+\beta_1\bar x+\bar\varepsilon, \qquad \bar\varepsilon=\frac{1}{n}\sum_{i=1}^n \varepsilon_i. \] Subtract the mean equation from the individual equation: \[ y_i-\bar y = \beta_1(x_i-\bar x)+(\varepsilon_i-\bar\varepsilon). \] Plug this into the expression of $\hat\beta_1$:
Now simplify the numerator: \[ \sum_{i=1}^n (x_i-\bar x)(\varepsilon_i-\bar\varepsilon) = \sum_{i=1}^n (x_i-\bar x)\varepsilon_i - \bar\varepsilon\sum_{i=1}^n (x_i-\bar x). \] But $\sum_{i=1}^n (x_i-\bar x)=0$, hence \[ \sum_{i=1}^n (x_i-\bar x)(\varepsilon_i-\bar\varepsilon) = \sum_{i=1}^n (x_i-\bar x)\varepsilon_i. \] Therefore the key decomposition is \[ \boxed{ \hat\beta_1 = \beta_1 + \frac{\sum_{i=1}^n (x_i-\bar x)\varepsilon_i} {\sum_{i=1}^n (x_i-\bar x)^2} }. \]
Take expectations conditional on $X=(x_1,\dots,x_n)$. Under the exogeneity assumption $E(\varepsilon_i\mid X)=0$, we obtain \[ E(\hat\beta_1\mid X) = \beta_1 + \frac{\sum_{i=1}^n (x_i-\bar x)E(\varepsilon_i\mid X)} {\sum_{i=1}^n (x_i-\bar x)^2} = \beta_1. \] Hence, by iterated expectations, \[ E(\hat\beta_1)=E\bigl(E(\hat\beta_1\mid X)\bigr)=\beta_1. \]
Recall the decomposition (simple regression with intercept): \[ \hat\beta_1 = \beta_1 + \frac{\sum_{i=1}^n (x_i-\bar x)\varepsilon_i}{\sum_{i=1}^n (x_i-\bar x)^2}. \] Define the shorthand \[ S_{xx}=\sum_{i=1}^n (x_i-\bar x)^2. \] Then \[ \hat\beta_1-\beta_1=\frac{1}{S_{xx}}\sum_{i=1}^n (x_i-\bar x)\varepsilon_i. \]
Compute the variance conditional on $X=(x_1,\dots,x_n)$: \[ \operatorname{Var}(\hat\beta_1\mid X)=\operatorname{Var}\!\left(\frac{1}{S_{xx}}\sum_{i=1}^n (x_i-\bar x)\varepsilon_i\;\middle|\;X\right) =\frac{1}{S_{xx}^2}\operatorname{Var}\!\left(\sum_{i=1}^n (x_i-\bar x)\varepsilon_i\;\middle|\;X\right). \]
Assume (i) homoskedasticity and (ii) no conditional correlation: \[ \operatorname{Var}(\varepsilon_i\mid X)=\sigma^2, \qquad Cov(\varepsilon_i,\varepsilon_j\mid X)=0\quad (i\neq j). \] Then
Therefore, \[ \boxed{ Var(\hat\beta_1\mid X)=\frac{\sigma^2}{S_{xx}} }. \]
By the law of total variance, \[ Var(\hat\beta_1)=E\!\left[Var(\hat\beta_1\mid X)\right]+Var\!\left(E(\hat\beta_1\mid X)\right). \] Under exogeneity, $E(\hat\beta_1\mid X)=\beta_1$, so the second term is zero and \[ Var(\hat\beta_1)=E\!\left[\frac{\sigma^2}{S_{xx}}\right]. \] For consistency, it is typically enough to work with the conditional variance.
Rewrite $S_{xx}$ as \[ S_{xx}=n\cdot \frac{1}{n}\sum_{i=1}^n (x_i-\bar x)^2. \] Under random sampling with $Var(x_i)=\sigma_x^2\in(0,\infty)$, we have the probability limit \[ \frac{1}{n}\sum_{i=1}^n (x_i-\bar x)^2 \xrightarrow{p} \sigma_x^2, \] so \[ S_{xx} \xrightarrow{p} n\sigma_x^2. \] Hence, \[ Var(\hat\beta_1\mid X)=\frac{\sigma^2}{S_{xx}} \approx \frac{\sigma^2}{n\sigma_x^2}, \] which implies \[ \boxed{ Var(\hat\beta_1\mid X)\xrightarrow{p} 0 \quad\text{and in particular}\quad Var(\hat\beta_1)=O\!\left(\frac{1}{n}\right)\to 0. } \]
Interpretation. As the sample grows, $S_{xx}$ increases proportionally to $n$ (provided $x$ has nonzero variance), so the sampling variability of $\hat\beta_1$ shrinks at rate $1/n$. This is the core variance argument behind the consistency of OLS slope estimates.
For a simple\footnote{Here, we omit the intercept to simplify the mathematics, but in fact you need to work with an intercept and the mathematics is just slightly longer to write down.} regression model: $y_i = \beta x_i +\epsilon_i$, $T$ the length of our data set, our Ordinary Least Square programme is: $$ \min_\beta \sum_i (y_i - \beta x_i)^2$$
With first order conditions, we can find the estimate of $\beta$ and we could prove that it is linear, unbiased, and efficient:
We can demonstrate that, if $\mbox{var}(\epsilon) = \sigma^2$: $$\mbox{var}(\beta) = \frac{\sigma^2}{\sum_i x_i^2}$$
The best estimate for the standard deviation of the residuals is: $$\hat{\sigma}^2 = \frac{\sum_i(\hat{y_i}-y_i)^2}{T-2}$$
And the standard deviation for the estimated parameters are\footnote{more details can be found \href{http://www.stat.cmu.edu/ hseltman/309/Book/chapter9.pdf}{here}} $$\hat{\sigma_{\beta_0}}^2 = \hat{\sigma}^2 \frac{\sum_i xi^2}{T \sum_i x_i^2 - \left(\sum_i x_i\right)^2}$$ $$\hat{\sigma_{\beta_1}}^2 = \hat{\sigma}^2 \frac{T}{T \sum_i x_i^2 - \left(\sum_i x_i\right)^2}$$
By decomposing the variance and with enough observations, it can be demonstrated that the total variance is decomposed by the variance of the model and the error terms: $$\mbox{var}(y_i) = \mbox{var}(\hat{y_i}) + \mbox{var}(\hat{\epsilon_i})$$ we then define $R^2$ as the ratio between the variance of the model and the total variance: $$R^2 = \frac{\sum_i (\hat{y_i}-\bar{y})^2}{\sum_i (y_i - \bar{y})^2} = \frac{\sum_i (\hat{y_i}-\bar{y})^2}{\sum_i (\hat{y_i} - \bar{y})^2 + \sum_i (y_i - \hat{y_i})^2}$$
$R^2 \in [0,1]$, the closest to $1$, the more explanatory our model (if our econometric methods were correct\footnote{for example, if the series have not a zero mean and we nonetheless force the intercept to be null, then the $R^2$ will seems larger than what the explanatory power of the model}), meaning the ability of our model to explain the movements of the dependent variable.
The t statistic for a coefficient is defined as: $$t = \frac{\hat{\beta}}{\hat{\sigma}_{\hat{\beta}}}$$
This statistics follow a Student law with $N-k$ degree of freedom ($N$ being the number of observations and $k$ the number of estimated parameters).
The null hypothesis $H_0$ is : $\beta = 0$, we can use a threshold of $5\%$, as it is two-sided, we search for $t^*$ such that: $$ \mbox{Prob}(X < t^*) = 0.025$$ If $\vert t \vert > t^*$, then we can reject $H_0$, again, if our econometric methods were correct, $x$ has a significant influence over $y$.
\quad
For a model with an intercept, the results become:
$$\hat{\beta_0} = \bar{y} - \hat{\beta_1}\bar{x}$$
And we can show that: $$R^2 = corr(yi,xi)^2$$
\stepcounter{question}
We might want to link the weight of a person with her height:
As we will describe in more detail section (ref), we do not seem to validate the homoscedasticity hypothesis, we take the log of weight so our data set better fit this hypothesis:
\stepcounter{question}
We assume that $\log(w)$ and $h$ follow given Data Generating Processes and might want to fit a linear regression:
From the literature on the subject, we know we are missing explanatory variables (parents' height and weight at a minimum, but we do not have that information).
Linear regressions are powerful tools, but they often stumble when the relationship between variables isn't strictly linear or when the data is highly skewed. One of the strategies to mitigate such issues is the log-linear regression.
The logarithm is a monotonic and concave function. This concavity property aids in reducing the skewness in datasets (for further details, refer to Section (ref)). By transforming the dependent variable with a logarithm, we can often achieve a more linear relationship with the predictor and ensure that the assumptions of linear regression are better met. This improved linearity often translates into more efficient and robust regression estimates.
For our study of the relationship between weight and height, the log-linear regression is particularly useful. The relationship can be expressed as:
Where:
Illustration: We hypothesize that the logarithm of an individual's weight, $\log(w_i)$, is linearly related to their height, $h_i$. Our model is: \[ \log(w_i) = \alpha + \beta h_i + \epsilon_i \]
This model suggests that the percentage change in weight is linearly related to changes in height.
cohn2022count in their Takeaway 2, demonstrate that bias due to heteroskedastic model error can cause log-linear regression coefficients to have the wrong sign. For this, they suppose that $y=\exp\left( \beta x \right) \eta$ where $x$ is normally distributed and $\eta$ is log-normally distributed with a standard deviation depending on the realization of $x$ (this is the modeling of heteroskedasticity), $\sigma_\eta(x)=\exp\left( \delta x \right)$ with $\delta$ a constant and they show that if $\delta > 2 \beta$ then this bias can cause the sign flip.
Note on Model Completeness: From prior research, we understand that crucial explanatory variables, such as parents' height and weight, influence an individual's weight. However, in this model, we don't account for these variables due to a lack of data. Such omission may lead to omitted variable bias.
Given the regression output:
- Intercept ($\alpha$): 4.2587. This is the expected logarithm of weight for an individual with a height of 0 (which is not meaningful in real-world context but provides a baseline for the model).
- Coefficient of $h_i$ ($\beta$): $0.0013$. Because $\tfrac{d \log w_i}{dx} = \tfrac{1}{w_i}\tfrac{dw_i}{dh_i}$, $\beta$ measures the proportional change in $y$: a one-unit increase in $h_i$ increases $w_i$ by about $100\times \beta = 0.13\%$, holding other factors constant.
- R-squared: 0.197. This means only 19.7% of the variability in the logarithm of weight is explained by height based on this model.
Given the statistically significant p-values for both intercept and slope, we infer that the relationship between height and the logarithm of weight is statistically significant. However, with an R-squared value of 0.197, height explains a relatively small fraction of the variation in the logarithm of weight.
When a data set has many 0s (e.g. $69\%$ fo Compustat firms are granted zero patents in a given year). To take those 0s into account, one might want to do a log1plus regression, because then 0s remain 0s with the log transformation.
cohn2022count demonstrates that the log1plus estimates is linked to the log-linear estimates by $$\lambda = \frac{E\left[ y\vert x \right]}{1+E\left[ y\vert x \right]} \beta$$ When $E\left[ y\vert x \right]$ is large then the two estimates are almost equal, but the authors note that when $E\left[ y\vert x \right]$ is large then it is unlikely that many observations are 0s.
For illustration only, here we perform a log1plus regression and as $E\left[ y\vert x \right]$ is large, then the two specification, log1plus and log-linear are similar and log1plus here doesn't add any value as we do not have any real $0$ in our data set.
While linear and log-linear regressions are commonly used for continuous data, count data or non-negative integer values often demand different models. One such model is the Poisson regression, typically used for count data.
In this context, citing cohn2022count, we introduce the Poisson regression mainly for illustrative purposes. It's worth noting that while our data of weight isn't typical 'count' data, the Poisson regression can still provide valuable insights.
The specification of the Poisson regression model is:
Where:
It's crucial to interpret the coefficients of Poisson regression in terms of percentage change, given its multiplicative nature.
bellego2022dealing provide some critics\footnote{(1) Existence of a solution is not guaranteed; (2) They can exhibit convergence issues; (3) They can be very imprecise if the log-scale error is heavy-tailed; (4) They can be difficult to estimate with many fixed-effects; (5) instrumental variables require stronger assumptions and may dramatically increase computational complexity; and (6) they suffer from an incidental parameter problem in the standard errors when including two-way fixed effects.} of the Poisson approach and suggest an alternative iOLS that goes back to log (one plus)-linear regression.
We can also use machine learning techniques, but these will be oriented toward answering the question: given a height $h_i$ what is the predicted weight $w_i$, but the model being nonparametric, we cannot easily derive a sensitivity.
As an illustration, we fit a Random Forest regression model on the sample, we might want to restrict our fit not to overfit:
knowing its limitation, the linear model allows an interpretation of the impact of height on weight.
Stationarity, Dickey-Fuller test (original and augmented), auto-regressions\footnote{stationarity_adf_ar.py}.
Before introducing the theory of stationarity and the estimation of AR, MA, ARMA, ARIMA, SARIMA, ARIMA--GARCH, and SARIMA--GARCH models, we address a set of frequently asked questions.
\noindentFormal decomposition of a SARIMA--GARCH model.
Let $\{y_t\}$ denote a univariate time series. A SARIMA$(p,d,q)\times(P,D,Q)_s$--GARCH$(r,s)$ model can be written as \[ \Phi_P(L^s)\,\phi_p(L)\,(1-L)^d(1-L^s)^D y_t = \Theta_Q(L^s)\,\theta_q(L)\,\varepsilon_t, \qquad \varepsilon_t = \sigma_t z_t, \] where $z_t \sim \text{i.i.d. }(0,1)$ and $\sigma_t^2$ follows a GARCH process.
\paragraph{Definition (White noise).} A stochastic process $\{\varepsilon_t\}_{t\in\mathbb{Z}}$ is said to be a white-noise process with mean zero and variance $\sigma_\varepsilon^2$, denoted \[ \varepsilon_t \sim \text{WN}(0,\sigma_\varepsilon^2), \] if it satisfies:
\paragraph{Definition (Moving Average process).} A time series $\{y_t\}_{t\in\mathbb{Z}}$ is said to follow a moving-average process of order $q$, denoted MA($q$), if it can be written as \[ \boxed{\; y_t = \mu + \varepsilon_t + \sum_{j=1}^{q} \theta_j\,\varepsilon_{t-j}, \qquad \varepsilon_t \sim \text{WN}(0,\sigma^2). \;} \]
Before regressing (or for AR), we need to go step by step on the integration order and the cointegration\footnote{For the cointegration test, the null hypothesis is no cointegration} of the variables.
A time series \(\{y_t\}\) is said to be weakly stationary (or second-order stationary) if it satisfies the following conditions for all time points \(t\) and for all lags \(s\):
These conditions ensure that the time series has stable moments (mean and variance) across different time points and does not exhibit trends or seasonality in its first and second moments.
It can be shown\footnote{Using the law of large numbers and assuming that \(y_t\) is stationary.} that the following empirical moments are consistent estimators of the above moments when \(T \rightarrow \infty\):
The most basic weakly stationary processes are white noises and moving averages.
A time series \(\{y_t\}\) is said to be strictly stationary\footnote{In financial time series, strict stationarity is rare but can be taken as an assumption and then discussed with robustness checks.} if the joint distribution of \((y_t, y_{t+1}, \ldots, y_{t+k})\) is identical to \((y_{t+h}, y_{t+h+1}, \ldots, y_{t+h+k})\) for all time points \(t\) and for all lags \(h\) and \(k\), that is,
In other words, the distribution of the series does not change over time when shifted along the time axis, making the series consistent in distributional properties at all time points.
If \(\Delta\) is the difference operator defined as \(\Delta y_t = y_t - y_{t-1}\), \(\Delta^2 y_t = \Delta(\Delta y_t)\), etc., we say that \(y_t\) is integrated of order d or \(I(d)\) if \(y_t\), \(\Delta y_t\), \(\ldots\), \(\Delta^{d-1} y_t\) are non-stationary and \(\Delta^{d} y_t\) is stationary\footnote{The concept is detailed in Engle:1987aa. Definitions vary; for example, alexanderII defines integrated of order \(d\) as \(\Delta^{d} y_t\) having a stationary ARMA representation. Keep in mind that integration orders greater than two are not very relevant for economic applications.}.
The Wold Decomposition Theorem plays a pivotal role in the theory of stationary time series. This theorem essentially describes any stationary time series as a combination of two orthogonal components: a deterministic component and a stochastic component.
Let \(\{y_t\}\) be a zero-mean weakly stationary time series. According to the Wold Decomposition Theorem, \(\{y_t\}\) can be decomposed as
where
The first component \(\sum_{j=0}^{\infty} b_j \varepsilon_{t-j}\) is known as the innovations or moving-average representation of \(\{y_t\}\). The deterministic component \(\eta_t\) is typically of less interest and is often taken as zero, especially when the focus is on the stochastic properties of the time series.
More compact wording for part of this theorem:
$\epsilon \sim \mbox{WN}(0,\sigma^2)$ and $\sum_{i=0}^\infty \vert b_i\vert < \infty$.
\paragraph{What does it mean for an infinite sum to converge?}
An infinite sum \[ \sum_{j=0}^{\infty} a_j \] is said to converge if the sequence of partial sums \[ S_N \equiv \sum_{j=0}^{N} a_j \] approaches a finite limit as $N \to \infty$. Intuitively, this means that adding more and more terms eventually has a negligible effect on the total sum.
Deterministic intuition. Consider the geometric series \[ \sum_{j=0}^{\infty} \rho^j, \quad |\rho|<1. \] Each additional term becomes smaller and smaller, and the total sum converges to $1/(1-\rho)$. In contrast, if $|\rho|\ge 1$, the terms do not decay fast enough and the sum diverges.
Impulse Response Function aims at answering how much our process ${y_t}$ is affected at time $t$ by an innovation (a shock) that happened $h$ periods before $t$.
The Impulse Response Function traces the effect of a one-unit shock to \(\varepsilon_t\) on the time series \(y_t\) at various time horizons. Mathematically, the \(h\)-step-ahead impulse response is given by \[ IRF(h) = E\left( y_{t+h} \mid \varepsilon_t = 1, \varepsilon_{t-1} = 0, \varepsilon_{t-2} = 0, \ldots \right), \] where \(h\) represents the horizon of interest.
Graphically, the IRF represents the response of \(y_t\) to a one-unit impulse in \(\varepsilon_t\), holding all other past \(\varepsilon_{t-i}\) constant at zero, illustrating how an external shock or innovation permeates and affects the time series over subsequent periods.
From the Wold theorem, for a stationary process the innovation at time $t-h$ is $\epsilon_{t-h}$ and the Impulse Response Function if $\frac{\partial y_t}{\partial \epsilon_{t-h}}$ which is thus $b_h$.
Consider a linear time series model (like an ARMA model): \[ y_t = \phi(L) y_{t-1} + \theta(L) \varepsilon_t, \] where \(y_t\) is the time series, \(\phi(L)\) and \(\theta(L)\) represent polynomial lag operators, \(L\) is the lag operator, and \(\varepsilon_t\) represents white noise.
We consider a stationary AR(1):
with $0<\vert \theta \vert < 1$ and $\epsilon \sim \mathcal{N}(0,1)$. We can rewrite: $$y_t = \sum_{i=1}^{t-1}\theta^i \epsilon_{t-i} + \theta^t y_0$$ we have:
We say that an AR(1) is asymptotically stationary, and usually just say stationary.
Another way to look at it, is that if this process is to be weakly stationary, then, we need:
which has solution if $0<\vert \theta \vert < 1$, as a variance is necessarily positive. We can rewrite:
If $\theta = 1$ either we have no solution, either $\sigma = 0$, which would then mean that the error term is a constant $\forall t, \epsilon_t = 0$ and then $\forall t, y_t = y_0$ which is stationary but presents no interest for study.
From our initial exploration of the AR(1) model, we can derive certain insights, especially when applied to financial data like asset prices. First, to ensure that $E\left(y_t\right) \rightarrow 0$, we should refrain from using the price levels directly. If we were to use them, it might erroneously suggest a model where all prices converge to zero, implying a rather grim scenario where every company defaults. Instead of price levels, it's more practical and meaningful to work with price returns, thereby de-trending the data to respect the condition $E\left(y_t\right) \rightarrow 0$.
For the constant volatility assumption inherent in the AR(1) model, let's consider the CAC 40 index time series introduced in Section (ref). Upon visual inspection, one can observe periods of 'volatility clustering', which implies that volatility isn't constant through time. This observation suggests that the basic AR(1) model might not be sufficient to capture the underlying dynamics of such financial time series. To address this, more advanced models that account for evolving volatility, such as the ARCH (Autoregressive Conditional Heteroskedasticity) model, become necessary.
Equation (ref) can be written using the lag operator $L$ as \[ (1-\theta L)y_t = \varepsilon_t. \]
The root of the autoregressive polynomial $(1-\theta L)$ is $L=1/\theta$. If $|\theta|<1$, this root lies outside the unit circle, and the lag polynomial is invertible.\footnote{The polynomial $\phi(L)=1-\theta L$ is called the autoregressive polynomial of $\{y_t\}$.}
In that case, \[ (1-\theta L)^{-1} = \sum_{i=0}^{\infty} \theta^i L^i, \] and equation (ref) admits the moving-average representation \[ y_t = \sum_{i=0}^{\infty} \theta^i \varepsilon_{t-i}. \]
Since $\sum_{i=0}^{\infty} |\theta|^i < \infty$, the process $\{y_t\}$ has finite variance and is weakly stationary.
This generalizes to an AR($p$) process written in lag-polynomial form as \[ \Phi(L) y_t = \varepsilon_t, \qquad \Phi(L) \equiv 1-\phi_1 L-\cdots-\phi_p L^p . \] Let $\{\lambda_k\}_{k=1}^p$ denote the roots of the characteristic equation $\Phi(z)=0$. If \[ |\lambda_k|>1 \quad \text{for all } k=1,\ldots,p, \] then $\Phi(L)$ can be factorized as \[ \Phi(L)=\prod_{k=1}^p (1-\lambda_k L), \] and therefore \[ \Phi(L)^{-1}=\prod_{k=1}^p (1-\lambda_k L)^{-1}. \] Since $|\lambda_k|>1$ implies $|\lambda_k^{-1}|<1$, each factor admits a geometric-series expansion, \[ (1-\lambda_k L)^{-1}=\sum_{j=0}^{\infty} (\lambda_k L)^j, \] so the product is a well-defined power series in $L$, \[ \Phi(L)^{-1}=\sum_{i=0}^{\infty} a_i L^i. \] Moreover, the coefficients $\{a_i\}$ decay geometrically, so in particular \[ \sum_{i=0}^{\infty} |a_i|<\infty, \] and the AR($p$) process admits the MA($\infty$) representation $y_t=\sum_{i=0}^{\infty} a_i \varepsilon_{t-i}$, implying weak stationarity.
If we follow section (ref) but now set $\theta = 1$, $y_t = y_{t-1} + \epsilon_t$ which is a Random Walk and chose $y_0=100$. We can rewrite the process with $L$ the lag operator, $(1 - L) y_t = \epsilon_t$, the root of the equation $1 - x$ is 1, hence we say that this process has a unit root.
We can rewrite $y_t = \sum_{i=1}^{t-1}\epsilon_{t-i} + y_0$ and find that $Var\left(y_t\right) = t \sigma^2 \rightarrow \infty$ which is against our definition of weak stationarity.
\quad
\stepcounter{question}
\quad
In simple terms, a time series has a unit root, or follows an integrated process, if shocks have permanent effects on the level, but not on the rate of change in the variable. Indeed, if at a time $t_{\mbox{choc}}$ the time series experience a shock $\Delta = 50$, then the future level of $y_t$ will be impacted ($\forall t \geq t_{\mbox{choc}}$) but not the rate of change $\frac{y_{t+1} - y_t}{y_t}$.
For illustration, we apply a shock $\Delta$ at mid-period to our time series:
We can follow the unit root test introduced by David Dickey and Wayne Fuller in 1979.
In its simplest form, the Dickey-Fuller test the null hypothesis $H_0$: $\rho-1=0$ with a Student test over the equation, with no autoregressive process:
Which can be re-written:
The statistics $\frac{\hat{\gamma}}{\mbox{Standard Error}\left(\hat{\gamma}\right)}$ can be tested again the threshold of the standard t-distribution, this is discussed and confirm in MacKinnon2010 for the simple version of the test.
To find the critical value, you cant import the t-distribution function from the stats package and define a probability threshold and degree of freedom\footnote{for $268$ degree of freedom and probability $95\%$ the critical value is $1.65$}.
This test is not powerful enough to have full confidence as there is likely to be autoregressive component in time series: use an augmented test.
The Augmented Dickey-Fuller test can be used to test whether a series is stationary, and take autoregressive effects into account. $H_0$: $\gamma = 0$.
The formula being used for the augmented Dickey-Fuller test is of the form:
MacKinnon2010 discuss the number of lags, selected via information criteria methods, and discuss the test statistics for the "no-constant" ('n'), "no-trend" ('c'), and "with-trend" ('ct') version of this test with reference to the table to be used.
\quad
With the python function adfuller in statsmodels, you can chose the parameter of the regression:
\quad
The output of the test are:
\quad with the null hypothesis:
$H_0$: the time series has a unit root.
\quad
To apply this test, you need to have a "testing strategy" as suggested in Elder:2001aa which "exploits known information about the growth character of the variable under investigation". So:
Additional data, you can refer to the book Lutkepohl2004 and the detailed answer posted \href{https://stats.stackexchange.com/questions/44647/which-dickey-fuller-test-for-a-time-series-modelled-with-an-intercept-drift-and}{here}.
Elder:2001aa stress the importance of the correct model as
\quad
The monthly variables of unemployment, population and GDP are integrated of order 1.
If no obvious conclusion can be taken from observations of the data set, then we can resort to a thorough procedure for the augmented Dickey-Fuller test, with three competing models:
For each model, the lag order ($p$, $q$ or $r$) can be selected by t-values of information criteria and one should check that with the chosen lag order the residuals have no serial correlation.
Then one starts following Dolado:1990aa with the model (ref) and tests for the null of $\beta_1$ if the null is rejected, the augmented Dickey-Fuller test can be performed with this model. If the null is not rejected, then one move on to model (ref). Note that for the null tests we must use the asymptotic distributions tabulated in Dickey:1981aa.
With the model (ref), one tests for the null of $\beta_0$, if the null is rejected, the aumented Dickey-Fuller test can be performed with this model. If the null is not rejected, then one move on to model (ref) and perform the augmented Dickey-Fuller test with this model.
We can define two time series as in Granger:1974aa: $$y_t = y_{t-1} + \epsilon_t$$ $$x_t= x_{t-1} + \nu_t$$ with $\epsilon_t$ and $\nu_t$ two independent white noise from a normal law and we set $y_0=0$ and $x_0=0$. When regressing: $$y_t = \alpha + \beta x_t + \gamma_t$$ one would expect the $\beta$ not to be significantly different than $0$ and the $R^2$ of the regression to be low. This is not the case as both series present a unit root.
The OLS estimator for \( \beta \) is: \[ \hat{\beta}_{OLS} = (X'X)^{-1}X'y. \] Given the non-stationary nature of both series the quadratic mean grows with time (see also the formula section (ref) to get a sense of why this variance tends to infinity with a unit root): \[ E(y_t^2) \propto T, \quad E(x_t^2) \propto T. \] Thus, as \( T \rightarrow \infty \), the variances of both series grow without bound. The covariance between the two series will also grow with \( T \) due to their shared trending behavior, affecting the \( X'y \) term. Given the nature of the OLS estimator: \[ \hat{\beta}_{OLS} = (X'X)^{-1}X'y, \] the matrices \( X'X \) and \( X'y \) grow to infinity due to the non-stationary nature of \( x_t \) and \( y_t \). This growth can lead to non-convergent behavior of \( \hat{\beta}_{OLS} \) even if the true \( \beta \) is zero (also get back to equation (ref) to get a sense of this result).
Moreover, the variance of the OLS estimator (given by the diagonal elements of the matrix \( \hat{\sigma}^2 (X'X)^{-1} \), where \( \hat{\sigma}^2 \) is the variance of the error term) can also be impacted by the non-stationarity of the series.
For the detailed of how to compute both the estimator and its variance, cf section (ref), I implement it directly in the code to illustrate how the numerator and/or the denominator can 'explode'.
phillips1986understanding demonstrates that he OLS estimator does not converge in probability as the sample size increases, the t- and F-test statistics do not have well-defined asymptotic distributions and the Durbin-Watson statistic converges to zero. The reason is that, with $y_t$ and $x_t$ being I(1) variables, the error term $\epsilon_t$ will also be a nonstationary I(1) variable.
\quad
We illustrate this with an example of two draws and their $R^2$, coefficients, t-values and Durbin-Watson statistics:
The $R^2$ seems "attractive"\footnote{and the time series have a correlation of $43\%$}, both the intercept and $\beta$ seems to have high t-value which means that they seems significantly different than $0$, but there is a low Durbin-Watson statistics which suggest that the residuals of the regressions are strongly autocorrelated. These are the symptoms of two I(1) series not cointegrated being regressed on one another, a spurious regression where two independent time series seems related. So in fact the true value of $\beta$ is zero but as demonstrated by Philips (1986) the estimate of $\beta$ doesn't converge in probability to zero (estimates of the regression coefficients are inefficient). Also the usual significance tests on the coefficients are invalid. This condition (high $R^2$ and low Durbin-Watson statistics) is not a necessary condition for a regression to be spurious as demonstrated in Granger:1974aa.
Autocorrelation, specifically serial correlation, refers to the correlation of a variable with itself through different time lags. This phenomenon violates the ordinary least squares (OLS) assumption that the error terms are uncorrelated, leading to inefficiency in estimators, and potentially distorting statistical inference.
The Durbin-Watson (DW) test is a formal statistical test of whether autocorrelation is present in the residuals of a regression analysis. It is commonly used in time-series data to detect the presence of serial correlation.
The Durbin-Watson statistic is calculated using the following formula: \[ d = \frac{\sum_{t=2}^{T} (\hat{\epsilon}_t - \hat{\epsilon}_{t-1})^2}{\sum_{t=1}^{T} \hat{\epsilon}_t^2} \] where:
The test statistic \(d\) usually lies between 0 and 4, where:
Specific critical values for the DW statistic can be found in statistical tables, which can be used to more formally accept or reject the null hypothesis of no autocorrelation.
Spurious regression often occurs when working with time-series data. It refers to the scenario where two or more time series variables are correlated due to shared trends over time, rather than because there is a genuine relationship between them.
A classic symptom of spurious regression is a very high \(R^2\) value along with significant t-ratios, despite no theoretical or logical connection between the variables. Additionally, residuals from spurious regressions often exhibit patterns, violating the assumption of no autocorrelation and potentially leading to the rejection of the null hypothesis in the Durbin-Watson test even in the absence of a meaningful relationship between the variables.
Identifying and rectifying spurious regressions typically involves:
Consider the AR(2)
where $\varepsilon_t$ is white noise (hence stationary).
\paragraph{Lag-polynomial and characteristic roots.} Using the lag operator $L$, (ref) becomes
Define the lag roots $\lambda_1,\lambda_2$ as the solutions to
Equivalently, \[ -\beta_2 \lambda^2-\beta_1 \lambda+1=0 \quad\Longrightarrow\quad \lambda_{1,2}=\frac{-\beta_1\pm\sqrt{\beta_1^2+4\beta_2}}{2\beta_2}\qquad (\beta_2\neq 0). \] Then the factorization is
\paragraph{Stationarity condition (and the role of $|\lambda|$).} A necessary and sufficient condition for (covariance) stationarity of the AR(2) part is
i.e. all lag roots lie outside the unit circle. Under (ref), the AR admits a convergent MA($\infty$) representation (mean-square convergence).
\paragraph{Companion matrix / eigenvalues.} Let $x_t=(y_t,y_{t-1})'$. Then \[ x_t=
+ A x_{t-1} +
, \qquad A=
. \] The eigenvalues $r_1,r_2$ of $A$ solve the characteristic equation
Stationarity is equivalently
The two parametrizations are consistent because $\lambda_i=1/r_i$: “$|\lambda_i|>1$” (lag roots) $\iff$ “$|r_i|<1$” (eigenvalues).
\paragraph{Unit root condition.} The process has a unit root at $L=1$ if and only if
since $L=1$ is a root of $1-\beta_1 L-\beta_2 L^2$ exactly under (ref). In that case the AR(2) is not stationary (it is $I(1)$ unless there is a second unit root).
\quad
\stepcounter{question}
\paragraph{Example $\beta_1=1.6$, $\beta_2=-0.6$.} Here $\beta_1+\beta_2=1$, hence there is a unit root. Indeed the eigenvalues of $A$ solve $r^2-1.6r+0.6=0$, giving $r_1=1$ and $r_2=0.6$, hence $\lambda_1=1$ and $\lambda_2=1/0.6\approx 1.67$. Because $|\lambda_1|=1$, condition (ref) fails: the level $y_t$ is not stationary.
\paragraph{First difference under a single unit root.} If (ref) holds but the other root satisfies $|\lambda_2|>1$ (equivalently $|r_2|<1$), then \[ (1-L)(1-r_2 L)y_t=\alpha+\varepsilon_t \quad\Longrightarrow\quad (1-r_2 L)\Delta y_t=\alpha+\varepsilon_t. \] Thus $\Delta y_t$ is AR(1) (with drift $\alpha$) and is stationary around a constant mean provided $|r_2|<1$.
\stepcounter{question}
Log transformations are often used in finance to stabilize variance and interpret relationships in terms of percentage changes. Key points include:
Log transformations are especially useful for series with multiplicative relationships or exponential growth but should be avoided if the data contains zeros or negative values.
If we believe\footnote{ARp_French_population.py} that that the time series of the French population is I(1) and there is no seasonality, we would want to fit an AR(p) process on the French population growth rate.
A first naive approach is to fit an AR(1), we de-trend the growth rate, $\tilde{y}_ t = y_t-\bar{y}$, with $\bar{y}$ the average growth rate observed over the period and then perform an ordinary least squares estimation of the linear regression: $$\tilde{y}_t = \theta \tilde{y}_{t-1} + \epsilon_t$$
As we are working with a finite sample of size $T$ and exogeneity can be safely assumed for shocks on past population growth rates but the exogeneity assumption doesn't hold for future population growth rate, then the OLS estimator $\hat{\theta}$ of $\theta$ is biased, $E\left(\hat{\theta}\right) - \theta \neq 0$. The demonstration will not be provided here, but you could try to convince yourself with the Example 2 of \href{http://web.econ.ku.dk/metrics/Econometrics2_05_II/LectureNotes/regression.pdf}{this lecture notes}. To put it simply: $$\hat{\theta} = \frac{\sum_{t=1}^T y_t y_{t-1}}{\sum_{t=1}^T y_{t-1}^2}$$ is unbiased if $E\left[ \epsilon_t\vert y_0, \ldots, y_t, \ldots, y_T\right] = 0$ which for time series is not the case.
But it can be shown that in case $\vert \theta \vert < 1$, then the series is strictly stationary and ergodic and the the OLS estimator is consistent, meaning that as the size of the sample increases (as $T \to \infty$) the estimate tends to the quantity being estimated.
Log-likelihood estimate of an AR(1) if we assume that errors are normally distributed:
If we assume that $y_t-\theta y_{t-1} \sim \mathcal{N}\left( 0, \sigma^2 \right)$ and
For the first observation $y_0$, we have a specific estimation of its distribution: $y_0 \sim \mathcal{N}\left(0, \frac{\sigma^2}{1-\theta^2} \right)$ (remember that $Var(y_t) = \theta^2 Var(y_{t-1}) + \sigma^2$).
as seen in section (ref) on maximum-likelihood estimation, the log-likelihood is additive $$\log\left(L(\theta)\right) = -\frac{T-1}{2}\log(2\pi\sigma^2)-\frac{1}{2\sigma^2}\sum_{t=1}^{T-1}\left(y_t - \theta y_{t-1} \right)^2 -\frac{1}{2}\log(2\pi) -\frac{1}{2}\log\left( \frac{\sigma^2}{1-\theta^2} \right) - \frac{1}{2}\frac{y_0^2}{\frac{\sigma^2}{1-\theta^2}}$$
We implement the log-likelihood estimate for this AR(1) on the French population and compare this estimate with the parameter from the OLS.
We might want here to try manually to fit a AR(3) and look at the significance of the regression coefficients: $$\tilde{y}_t = \theta_1 \tilde{y}_{t-1} + \theta_2 \tilde{y}_{t-2} + \theta_3 \tilde{y}_{t-3} + \epsilon_t$$
We will see starting section (ref) how to select the order $p$ of the AR, also look at moving average (MA) component and how to fit the model to the data. Methods used are \href{https://math.unice.fr/ frapetti/CorsoP/Chapitre_4_IMEA_1.pdf}{maximum likelihood estimations and method of moments}.
In simple terms, the Partial Autocorrelation Function (PACF) of order $k$ can be estimated as $P_{k}$ in the following AR(k) estimations: $$\tilde{y}_t = P_{1} \tilde{y}_{t-1} + \epsilon_{1,t}$$ $$\tilde{y}_t = \theta_{1,2} \tilde{y}_{t-1} + P_{2} \tilde{y}_{t-2} + \epsilon_{2,t}$$ $$\tilde{y}_t = \theta_{1,3} \tilde{y}_{t-1} + \theta_2 \tilde{y}_{t-2} + P_{3} \tilde{y}_{t-3} + \epsilon_{3,t}$$ $$ \ldots $$
The significance level at $95\%$ for a series with $T$ observations is\footnote{which can be approximated directly with $\pm \frac{2}{\sqrt{T}}$} $\pm \Phi^{-1}\left(\frac{1+0.95}{2}\right) \frac{1}{\sqrt{T}}$ with $\Phi^{-1}$ the inverse of the cumulative distribution function of the normal distribution.
To conclude on this AR models, to chose the order $p$, you can introduce and minimize the Akaike's Information Criterion (AIC) which you compute at each order $p$: $$AIC_p = - \frac{2}{T} log\left(L(\theta,p)\right) + \frac{2}{T} p $$
or the Schwarz-Bayesian information criterion (BIC): $$BIC_p = \log\left(\hat{\sigma}^2(p)\right) + \frac{\log(T)}{T} p $$
with $\hat{\sigma}(p)$ the estimate of the innovation standard deviation from the maximum likelihood of the AR at order $p$. BIC tends to penalize more lag order $log(T)$ for large data set.
We usually strive for model parsimony (limiting the order $p$ in our example), the first reason is the Occam's razor principle where the least assumptions you make the more likely your explanation. The second is the risk of overfitting where if you have $T$ observations, you could reach $p \to T$ and we can mention in that sense Ledolter:1981aa that argued that each unnecessary parameter increases the variance of the prediction error by a factor of $\frac{\sigma^2}{T}$.
One might wonder which criteria to use between AIC and BIC, but we are not aware of one that outperform the other in every applications. BIC might lead to model that are too parsimonious, this is why as a practitioner we compare both. We discuss further lag order selection in section (ref).
We use the procedure to fit an ARMA (although restricting MA order $q$ to null for this exercise).
If our assumptions are correct, the residuals from the model should be a white noise. We apply a Ljung-Box test\footnote{see section (ref) for more detail on this test} removing the order $p$ of the AR to the degree of freedom of the chi-square distribution.
To test for serial correlation with a maximum lag $h$ in the residual time series returns, we compute the Ljung-Box q-statistic against a $\chi_h^2$ distribution as in LJUNG:1978aa:
with $\tau_k$ the sample autocorrelation at lag $k$.
For residuals of an AR model, we remove the order $p$ of the AR to the degree of freedom of the chi-square distribution.
The $H_0$ of the Ljung-Box test is: no serial correlation at lag max.
With first a visual inspection of the AR(3) model residuals, they don't seem to be iid:
With this test, we reject the hypothesis that residuals are independently distributed, but we conclude that residuals are rather serially correlated.
Our model is not satisfactory.
In fact, seasonality for birth rate in France is a well documented fact \href{https://www.ined.fr/fichier/s_rubrique/168/pop_soc474.fr.fr.pdf}{as INED showed}, as such the 23rd of September has $5\%$ additional births as a usual days, this day is exactly 256 days after the 31st of December.
The main question is: how to treat seasonality (in French population, in a price index)? And especially, should you first treat seasonality and then do the first difference to have a time series stationary. We know that a time series with seasonality cannot be stationary, so we would like to think that we first need to treat seasonality and then make the time series stationary. But in fact, my take is that you cannot do this sequentially if you want to do it properly from an econometric point of view. This is that same as for the question of the ARMA(p,q)-GARCH(r,s) section (ref). You have to fit with MLE p, q, r, and s at the same time and not do first an ARMA and then fit a GARCH on the residuals. This is why in fact I would advise to learn time series until you master SARIMA and then treat a univariate time series that exhibit seasonality with an SARIMA. My advice would then be:
If the original time series data you work with has trend\footnote{to deal with trend, a first difference or computing the rate of change should do} and/or seasonality, then this series is unlikely to be stationary (or the power of your ADF test can be questioned). You will want to deal with this seasonality and trend before you try to fit a model. But even before this, if your time series display increasing variability, you need to address this.
In this section we follow\footnote{seasonality_vansteenberghe.R}, the material from the books Econometric models and economic forecasts from Robert Pindyck and Daniel Rubinfeld published in 1998 and Analyse des series temporelles from Regis Bourbonnais and Michel Terraza published in 2010.
Starting with the artificial data from Analyse des series temporelles from Regis Bourbonnais and Michel Terraza published in 2010, we observe three problems: increasing variability, trend and seasonality.
We first need to tackle the increasing variability issue, with a Box-Cox transformation, the most famous of which is the $\log()$.
Then we deal with the monthly seasonality with a difference with 12 lag, this also takes care of the trend.
If the time series present some long term non-linear behaviour (e.g. logarithm, exponential, power), we can use a Box-Cox transformation to revert to a linear process. The transformation can be written:
The purpose of a Box-Cox transformation is to modify the data set so it is more normally distributed.
\stepcounter{question}
We focus on the coffee world prices with 240 observation for Brazil downloadable \href{http://www.economicswebinstitute.org/data/coffeeprices.zip}{here}.
We get the following Buys-Ballot table, it is not obvious to observe month or year effect:
\resizebox{\textwidth}{!}{
}
We have a data set with:
We can define the following "means": $$x_{i.}=\frac{1}{p}\sum_{j=1}^p x_{ij}$$ $$x_{.j}=\frac{1}{N}\sum_{i=1}^N x_{ij}$$ $$x_{..}=\frac{1}{N p} \sum_{j=1}^p \sum_{i=1}^N x_{ij}$$
The following sum of squares: $$S_p=N \sum_{j=1}^p (x_{.j}-x_{..})^2$$ $$S_A=p \sum_{i=1}^N (x_{i.}-x_{..})^2$$ $$S_r=N \sum_{j=1}^p (x_{ij}-x_{i.}-x_{.j}+x_{..})^2$$
The following variances: $$V_p=\frac{S_p}{p-1}$$ $$V_A=\frac{S_A}{N-1}$$ $$V_r=\frac{S_r}{(p-1)(N-1)}$$
\quad
We cant compare the Fisher statistics: if the empirical Fisher statistics $$F_c = \frac{V_p}{V_r}$$ is greater than the Fisher statistics $F_{(p-1) , (N-1)(p-1)}^{0.05}$ then we can reject $H_0$ and the series is seasonal. if the empirical Fisher statistics $$F_c = \frac{V_A}{V_r}$$ is greater than the Fisher statistics $F_{(N-1) , (N-1)(p-1)}^{0.05}$ then we can reject $H_0$ and the series has a trend.
We have monthly data, for each observation we create a dummy, which will be equal to $1$ when the month of the observation is that given month, $0$ otherwise. Testing for seasonality is testing for the month significance. None can be find here.
We use a data set from the \href{http://www.bom.gov.au/climate/change/index.shtml#tabs=Tracker&tracker=timeseries}{Autralian Governement - Bureau of Meteorology} which offers time series on monthly rainfall data in Australia:
When the seasonal component seems visually on the plot to depend on the level, we apply a multiplicative decomposition.
Multiplicative seasonal adjustment rely on the fact that each time series $y_t$ can be decomposed as:
with
The aim is to decompose the time series, isolate and combine $L_t \times C_t$ (estimate it) then we can deduct an estimate of $S_t \times I_t$.
If we have monthly data, we estimate:
Then:
Then $I$ has to be eliminated, we average each $\hat{S}_t \times \hat{I}_t $ element is average per month ($\tilde{z}_1$ for January, $\tilde{z}_2$ February, etc.). Those terms should sum to 12. If they do not, we define $\bar{z}_i=\tilde{z}_i \frac{12}{\sum_j \tilde{z}_j}$.
Finally, the seasonally adjusted series are:
For rain, as the seasonal component should not vary with time nor the level of rain, we tend to use the simplest additive model $Y_t = T_t + S_t + \epsilon_t$ that we can rewrite as $Y_t= L_t + S_t + C_t + I_t$ where $T_t= L_t + C_t$. To 'deseasonalize', we want to work with $Y_t - S_t$.
In order to separate the long-term secular trend from the long-term cyclical component, we can apply a Hodrick-Prescott filter on $T_t$, for monthly data we use $\lambda=129 600$.
R should allow the use of the United Sates Census Bureau seasonal adjustment program, \href{https://www.census.gov/srd/www/x13as/}{X-13ARIMA-SEATS}.
We use\footnote{SARIMA_vansteenberghe.R} the data from the astsa package.
The dataset, graciously provided by Dr. Roy Mendelssohn of the Pacific Environmental Fisheries Group (personal communication), is part of the astsa package and spans 453 months from 1950 to 1987. It includes:
The objective is to fit an AR(2) model that accounts for seasonality and compare its performance against models that disregard seasonality, like a simple ARMA(3,1):
The SARIMA model, a multiplicative extension of the traditional ARIMA, accounts for both non-seasonal and seasonal components in the data. Formally, a SARIMA model is denoted as ARIMA$(p,d,q) \times (P,D,Q)_s$, and its mathematical representation is:
To determine the optimal values for the parameters, it's common to inspect the autocorrelation function (ACF) and the partial autocorrelation function (PACF)\footnote{\href{https://www.otexts.org/fpp/8/9}{Rob J Hyndman and George Athanasopoulos}}. The seasonal components of an AR or MA model manifest as lags in the PACF and ACF:
SARIMA models, given their ability to account for seasonal fluctuations, offer a robust approach to time series forecasting. By applying SARIMA to fisheries data, we aim to glean insights that are both economically significant and environmentally informative.
\stepcounter{question}
Varves are melt of glacier deposit. With the package asta, we have varves data in New England for 634 years. We can observe from the original series that the thicker the deposited amount, the higher the volatility, a logarithmic transform can correct this nonstationarity in the variance.
\stepcounter{question}
For this section, we use the detailed Gauss-Newton method presented in the book ngai.
An MA(1) can be written as: $y_t = \epsilon_t - \theta \epsilon_{t-1}$. If $\theta < 1$, then we can write: $\epsilon_n = y_n + \theta y_{n-1} + \ldots + \theta^{n-1} y_1$. Note that in this case, we have $n$ observations of our process $y$. We can search $\theta$ that minimizes: $S(\theta) = \sum_{t=1}^n \epsilon_t^2$.
We can compute this iteratively as we know that: $$\epsilon_1 = y_1$$ $$\epsilon_2 = y_2 + \theta \epsilon_1$$ $$\ldots$$ $$\epsilon_n = y_n + \theta \epsilon_{n-1}$$
We find that our $\theta$ lies between $-0.8$ and $-0.6$:
\stepcounter{question}
As observed in the ACF of the varves time series, there is a long memory in the process. If we want to avoid first differencing our data, we can revert to a fractional differenced series: $-.5 < d < .5$
Using the fracdiff function and fitting it on the varve data, we find a value of $d=0.38$
\stepcounter{question}
In this section\footnote{pandas_exercise_part_2_seasonality_French_population.py}, we focus on analyzing the monthly population data of France obtained from INSEE (Institut National de la Statistique et des Études Économiques). The dataset, identified as "Démographie - Population au début du mois - France métropolitaine" (ID: 000436387), provides a valuable insight into the demographic changes over time. Our objective is to investigate the presence of seasonality in this time series and apply appropriate econometric techniques to model this phenomenon.
Preliminary Data Analysis: The initial step involves loading and preprocessing the data using Python's pandas library. We ensure that the time series is indexed by the date, facilitating easier analysis. A preliminary visual inspection is conducted to identify any noticeable trends or seasonal patterns.
Seasonality in Population Data: Given the nature of demographic data, it is reasonable to hypothesize the presence of seasonal effects. For instance, certain months might exhibit higher population growth rates due to factors like birth rates, migration patterns, etc.
There are four main methods to deal with monthly seasonality in growth rates:
To assess the presence of seasonal effects, we initially employ a regression model incorporating dummy variables for each month. This method allows us to isolate the impact of each month on the population growth rate.
To extend our analysis, we conduct a comparative regression analysis considering each month as a categorical variable. Here, September serves as the reference category.
Upon examining the results from the regression models, we observe that:
Stationarity and Unit Root Testing: Before proceeding with the SARIMA model, it is crucial to ensure that the time series is stationary. We employ the Augmented Dickey-Fuller (ADF) test under different specifications (with constant, with constant and trend, and without both) to test for unit roots. The presence of unit roots would imply non-stationarity, necessitating differencing or transformation of the series.
SARIMA Model Specification: The SARIMA model, which stands for Seasonal AutoRegressive Integrated Moving Average, is a natural choice for modeling time series data with seasonality. The model is specified as SARIMA(p, d, q)(P, D, Q, S), where:
Parameter Selection and Model Tuning: Selecting the optimal parameters for SARIMA is done through a systematic grid search approach. We iterate over a range of values for each parameter and select the combination that minimizes an information criterion, typically the Akaike Information Criterion (AIC). This process helps in identifying a model that balances goodness-of-fit with parsimony.
The optimal SARIMA model identified for the monthly population growth rate data is SARIMA(0, 1, 0)x(1, 0, 1, 12). The model can be expressed as follows:
In this equation, \( Y_t' \) denotes the first-order differenced time series at time \( t \). The model includes a seasonal autoregressive term (\( \alpha Y_{t-12}' \)) and a seasonal moving average term (\( \theta \epsilon_{t-12} \)), both of which operate on a 12-month lag, reflecting the annual seasonal pattern in the monthly data.
Diagnostic Tests for SARIMA Model Residuals
We start by forecasting future population growth rates using our fitted SARIMA model. The model generates predictions for the growth rate from March to October 2020. This forecast is based on the underlying patterns and seasonal cycles identified in the historical data.
We then plot these forecasted growth rates alongside the actual growth rates observed during the same period to visually assess the model's accuracy.
To understand the impact of these growth rates on actual population numbers, we convert the predicted growth rates into estimated population levels. This is done by first adjusting the forecasted growth rates and then cumulatively multiplying these rates with the last known population level:
Similar adjustments are applied to the forecast's confidence intervals to estimate the range of possible future population levels.
Finally, we plot the projected population levels against the actual population levels observed from March to October 2020. This comparison provides a clear visual representation of the model's performance in predicting not just the rate of population growth but also the actual population size, or can help us understand the impact of the Covid-19 crisis.
Inflation\footnote{vansteenberghe_inflation_persistence.py} is one of the key variables in monetary economics, both for policy design and for forecasting. A central question concerns its persistence, that is, the long-run effect of a temporary shock to inflation: how strongly and for how long does a 1% increase in inflation today affect future inflation rates, and how quickly does inflation return to its previous level, if at all?
In this exercise, we study inflation persistence using monthly data from the Eurostat Harmonised Index of Consumer Prices (HICP), excluding tobacco. Following Chapter 8.8 of verbeek and the empirical approach in pivetta2007persistence, inflation dynamics are modeled using autoregressive processes.
The objective is threefold.
\paragraph{Data source.} The inflation series used in this exercise comes from Eurostat and corresponds to the Harmonised Index of Consumer Prices (HICP) for the European Union (EU aggregate). We focus on the headline index (all items), expressed as an index with base year 2015.
\paragraph{Seasonal adjustment.} Monthly inflation indices exhibit strong and recurrent seasonal patterns (e.g.\ sales periods, energy price seasonality). To isolate medium- and long-run movements relevant for persistence analysis, we remove the seasonal component using a standard seasonal decomposition.
We rely on a multiplicative decomposition rather than an additive one. For price indices, seasonal fluctuations tend to scale with the level of the series: a given percentage seasonal effect corresponds to a larger absolute change when the index level is higher. The multiplicative specification is therefore more appropriate and yields a seasonally adjusted series obtained by dividing the observed index by its estimated seasonal component.
\paragraph{About X-13ARIMA-SEATS.} As a robustness check, one could use the X-13ARIMA-SEATS procedure, the reference seasonal adjustment method employed by statistical agencies. X-13 combines ARIMA-based pre-adjustment, outlier detection, trading-day corrections, and signal extraction. While more powerful and flexible, it requires an external installation and is unnecessarily complex for a first exposure to seasonal adjustment. For pedagogical purposes, the simple decomposition approach is sufficient and fully transparent.
\paragraph{Unit root testing: ADF vs.\ KPSS.} Following verbeek, we rely primarily on the Augmented Dickey--Fuller (ADF) test, whose null hypothesis is a unit root. For inflation series, ADF results are often marginal and sensitive to the number of lags, so sharp conclusions are rarely robust. The KPSS test, which reverses the null (stationarity), typically leads to similarly ambiguous conclusions due to its dependence on bandwidth choices.
\paragraph{KPSS test and interpretation.} To circumvent the low power of standard unit root tests, kwiatkowski1992testing propose a test in which (trend-)stationarity is the null hypothesis and the existence of a unit root is the alternative. This test, known as the KPSS test, decomposes the observed series as \[ y_t = d_t + r_t + \varepsilon_t, \] where $d_t$ is a deterministic component (intercept or intercept and linear trend), $r_t$ is a random walk component, and $\varepsilon_t$ is a stationary error term, which is not required to be white noise. Under the null hypothesis of (trend-)stationarity, the variance of the random walk innovation is zero, so that $r_t$ is degenerate.
The KPSS test is a Lagrange multiplier test. Its computation proceeds as follows. First, an auxiliary regression of $y_t$ on $d_t$ is estimated by OLS and the residuals $\hat e_t$ are saved. Next, partial sums $S_t = \sum_{s=1}^t \hat e_s$ are formed. The test statistic is then given by \[ \mathrm{KPSS} = \frac{1}{T^2} \sum_{t=1}^T \frac{S_t^2}{\hat\sigma^2}, \] where $\hat\sigma^2$ is an estimator of the long-run variance $\sigma^2 = \sum_{j=-\infty}^{\infty} \mathbb{E}(\varepsilon_t \varepsilon_{t-j})$.
Estimation of $\hat\sigma^2$ requires the choice of a kernel (weighting scheme) and a bandwidth. Common choices include Bartlett weights, as in kwiatkowski1992testing, and the quadratic spectral kernel andrews1991heteroskedasticity. In practice, the KPSS test is highly sensitive to these choices, particularly in finite samples, and tends to over-reject stationarity when persistence is high.
Because the ADF and KPSS tests place opposite hypotheses in the null and both suffer from low power and tuning-parameter sensitivity in persistent macroeconomic series, they are best interpreted as complementary diagnostic tools rather than decisive tests. Accordingly, in this course, inflation persistence is quantified directly using autoregressive models, through the sum of AR coefficients and the implied half-life of shocks.
\paragraph{AR models in Stock and Watson (2007).} In stock2007has, the central univariate benchmark for inflation forecasting is a finite-order autoregressive model, denoted AR(AIC). Inflation is modeled in first differences to allow for a unit root in levels, and the lag order is selected recursively by the Akaike Information Criterion (AIC). This AR(AIC) specification serves as the baseline against which all alternative models are evaluated.
\paragraph{Method and lag order.} Operationally, the AR(AIC) model is estimated as a direct multi-step regression of future average inflation on lagged changes in inflation, with the number of lags chosen by minimizing AIC at each forecast origin. This approach allows for flexible short-run dynamics while remaining fully univariate and transparent. In the postwar U.S.\ data, Stock and Watson show that low-order AR representations provided a good approximation during the Great Inflation, but their forecasting performance deteriorated after the mid-1980s.
verbeek works on a "simple" AR model, citing the literature (older)
\paragraph{AIC vs.\ BIC.} Both the Akaike Information Criterion (AIC) and the Bayesian Information Criterion (BIC) trade off model fit against parsimony. AIC tends to select richer models, while BIC applies a stronger penalty that increases with sample size and favors more parsimonious specifications.
\paragraph{Practice in this exercise.} Here, both AIC and BIC select an AR(6) specification, which strengthens confidence in the chosen lag order. When criteria diverge, AIC is usually favored for inflation dynamics, while BIC provides a useful robustness benchmark.
Following Chapter 8.8 in verbeek and pivetta2007persistence, measure the inflation persistence. Use Eurostat Harmonised Index of Consumer Prices (HICP) excluding tobacco.
\stepcounter{question}
An $ARMA(p, q)$ model can be transformed into an $MA(\infty)$ model to facilitate the computation of forecast errors. The $ARMA(p, q)$ model is written as:
\[ y_t = \sum_{j=0}^{\infty} \alpha_j \varepsilon_{t-j} \quad \text{with} \quad \alpha_0 = 1 \]
where \( \alpha_j \) are the coefficients to be determined from the model. The \( h \)-period-ahead forecast, in terms of the errors \( \varepsilon_t \), is given by:
\[ y_{T+h\vert T} = E\left[y_{T+h}\vert y_T, y_{T-1}, \dots \right] = \sum_{j=h}^{\infty} \alpha_j \varepsilon_{T+h-j} \]
The variance of the forecast error can be expressed as:
\[ E \left[(y_{T+h} - y_{T+h\vert T})^2\right] = \sigma^2 \sum_{j=0}^{h-1} \alpha_j^2 \]
This formula allows us to compute the variance of the forecast errors for any \( h \)-period-ahead forecast. The variance increases with the number of forecast steps, reflecting the increased uncertainty over time. Using the forecast error variance, we can construct 95% confidence intervals for the forecasts. For a one-step-ahead forecast, the interval is given by:
\[ y_{T+1|T} \pm 1.96 \sqrt{c_1^2} \]
where \( c_1^2 = \sigma^2 \sum_{j=0}^{0} \alpha_j^2 = \sigma^2 \). For \( h \)-period-ahead forecasts, the confidence interval is given by:
\[ y_{T+h|T} \pm 1.96 \sqrt{c_h^2} \]
where
\[ c_h^2 = \sigma^2 \sum_{j=0}^{h-1} \alpha_j^2 \]
The wider the confidence interval, the greater the uncertainty in the forecast due to the accumulation of forecast errors as \( h \) increases.
For an AR(1) model, we have \( \alpha_j = \theta^j \). The expected quadratic forecast errors are given by:
\[ c_1^2 = \sigma^2, \quad c_2^2 = \sigma^2 (1 + \theta^2), \quad c_3^2 = \sigma^2 (1 + \theta^2 + \theta^4), \dots \]
For \( h \to \infty \), the forecast error variance converges to:
\[ c_{\infty}^2 = \frac{\sigma^2}{1 - \theta^2} \]
This represents the unconditional variance of \( y_t \), showing how forecast uncertainty increases over time but stabilizes in the long run for stationary processes.
As an illustration of in-sample "precision" of our model, we plot the normal distribution of the residuals and the shaded area equals $95\%$, we plot similarly the boundaries on the time series plot forecast:
For a $5\%$ two-sided confidence interval we take a t distribution with n degrees of freedom (which we approximate by a normal distribution as $n >> 120$): $Z_{\alpha/2}=1.96$.
\stepcounter{question}
Following barnard1963, plot the data of airline passengers and follow the Box-Jenkins method and the adaptive forecasting method.
We fetch\footnote{code: ARIMA_vansteenberghe.py} CAC 40 index historical data from 1995-1-16 to 2017-1-16. We compute the daily return on business days only.
One of the simplest question we could be interested in when it comes to financial time series is: "what is its value most likely to be tomorrow?"
If we assume (wrongly) that the daily returns are iid and follow a Gaussian DGP, then we only need to look at the histogram of the observed daily returns and can follow the logic presented section (ref), we model the CAC 40 return as a white noise section (ref).
But if we call $y_t = \frac{P_t}{P_{t-1}}-1$ the daily returns, we can expect that the realization of $y_t$ will depend on the realization of $y_{t-k} , k\in\mathbb{N}$.
One would take the historical data available for this time series, in our example, we have 5 638 observations from 1995-1-16 to 2017-1-16. We could choose the parameters based on the first $60\%$ of our data set (train our model) and then verify with the remaining $40\%$ observation that our model performs well (test our model). We might also want to cross-validate it, \href{http://scikit-learn.org/stable/modules/cross_validation.html}{more on this topic can be found here}.
A white noise is a is a sequence of independent and identically distributed random variables $y_t$ with finite mean and variance. Most of the time we take the assumption that $y_t$ is normally distribution
We use the simple model:
with $\epsilon_t \sim \mathcal{N}(0,\sigma^2)$ and as we are working on modelling and forecasting the CAC 40 returns: $\mu$ is the mean of the daily returns.
As defined in section (ref), this series is by definition stationary and as expected, when we perform an augmented Dickey-Fuller test: the series comes out as stationary.
We modify our model with a trend:
with $\epsilon_t \sim \mathcal{N}(0,\sigma^2)$.
We take $\delta=0.1\%$.
We perform an augmented Dickey-Fuller test taking the trend into account: the series is stationary.
If we are interested in the significance of the trend component of this test, we might have to do it by hand. Let's consider:
We find that in this case, the trend component is significant.
\quad
We might want to remove the trend. We can start with a first difference.
\quad
We also study quadratic trend model.
If we think that the CAC 40 returns are following a random walk, then:
s.t. $$ \epsilon \sim \mathcal{N}(0,\sigma^{2})$$
\quad
Here it is just that we use the previous business day's return to predict today's return: $E(y_t) = y_{t-1}$
This brings us to some of the measurements of forecasting errors:
We follow section (ref) that introduced AR(1) models.
Before starting, intuitively, an AR(1) model would mean that $y_t$ and $y_{t-1}$ related, therefore we can test for serial correlation of our daily returns, e.g. with a Ljung-Box test introduced section (ref).
There are two reasons to ditch the CAC 40 returns for an AR(1) model (no serial correlation, our series looks more like a withe noise) and the AR(1) parameters are not significant. For a further study of ARCH effect in CAC 40 returns, see section (ref).
We introduce a stationary AR(1) on time series of better returns candidates and for stationary, we check that: $0<\vert \theta \vert < 1$:
with $\epsilon_t \sim \mathcal{N}(0,\sigma^2)$, thus by definition (or assumption in our model if we apply it to a real time series) we assume that $y_t$ and $\epsilon_t$ are uncorrelated.
We can search for the optimal value of $\phi$, once we re-ordered the AR(1), detrending the $y_t$ time series to get rid of the intercept in the AR(1):
If $\phi<1$ we can iterate backward :
If we note $x_t=y_t-\mu$, $E_t(x_t) = \phi x_{t-1}$. That is why AR(1) is usually a benchmark model against which one assess the relevance of a more complicated model. If a complex model cannot beat an AR(1) forecast, then we might just stick to an AR(1) model to base our predictions on.
Hence the Autocorrelation Function (ACF) of an AR(1) is: $\rho(h)=\frac{\gamma(h)}{\gamma(0)}=\phi^h$, $h\geq 0$
If $\vert \phi \vert > 1$ the AR(1) process is explosive.
\stepcounter{question}
The definition of an AR(p) model is simply:
MA(1) model:
with $$E(y_t)=\mu$$ $$ \epsilon \sim \mathcal{N}(\mu,\sigma^{2})$$
If we note: $\gamma(h) = cov(y_{t+h},y_t)$
Hence the ACF of an MA(1) is:
We can demonstrate that there is a problem of non-uniqueness of a MA(1) model (take $\theta$ or $\frac{1}{\theta}$ and an adjusted value of $\sigma$ for example).
The definition of an MA model of order $q$ is simply:
The autocorrelation function (ACF) of a stationary time series is: $$\rho_k = corr(y_t,y_{t+k})$$
We call it a partial autocorrelation function (PACF) of a stationary process when we filter out the effect of the random variables $y_{t+1}, \ldots , y_{t+k-1}$. The PACF is used to identify the order ($p$) of an autoregressive model, i.e. AR($p$).
and
with $ \hat{y}_{t+k+1}$ (respectively $ \hat{y}_{t+1}$) being the orthogonal projection of $y_{t+k+1}$ (respectively $y_{t+1}$) on $y_{t+1}, \ldots , y_{t+k}$.
We can use the ACF and the PACF as visual tools to identify respectively the order of the Moving Average and the Auto Regressive one:
A process $y_t$ is an ARMA(p,q) if it is stationary and:
with $\phi_p\neq 0$ and $\theta_q \neq 0$
We can follow the Box-Jenkins methodology:
\quad
For our time series, we reject $H_0$ for the Breusch-Godfrey test and residuals seems to be autocorrelated, our model is not acceptable and we would need to search for a better model.
\quad
Topics to further explore: parameters redundancy, causality, invertibility
cairns2000discussion:
For the orders selection, section (ref) also provide some guidelines (in the VAR context).
Most time series can be decomposed in a nonstationary trend component and a zero-mean stationary composant: $$y_t=\mu_t +x_t$$ such that $$\mu_t=\beta_0 +\beta_1 t$$
$x_t$ is stationary thus $\Delta y_t = \beta_1 + \Delta x_t$ is stationary.
\quad
A process $x_t$ is said to be ARIMA(p,d,q) if $\Delta^d x_t$ is ARMA(p,q).
We follow\footnote{20190625_Box_Jenkins_methodo.R} the Box-Jenkins methodology to fit an ARIMA model as introdued in box1970.
Once we identified the orders p,d and q of the ARIMA, we estimate the parameters, e.g. using AIC criterion and we can simulate our model. Again, here we are not looking for a perfect reproduction of the historical time series, rather we want to check that our model can replicate the behavior of the original time series.
Once estimated, we check the error terms using the Ljung-Box\footnote{see section (ref) for more detail on this test} test to test the null hypothesis: $H_0$: the error terms are independent.
\stepcounter{question}
We want to forecast the the returns and for this explore a DGP that would be 'behind' the observations. If the observations were iid, then $E\left(y_t\vert y_{t-1} \right) = E\left(y_t \right) = \bar{y}$. Here this is not the case and a simple AR(1) model beat a forecast only based on the mean. We next need to explore the residuals of this AR(1) model and check if they are iid with a 'well behaved' DGP.
If we observe the residuals $\epsilon_t$ of an ARMA(1,1) model of the daily returns, we observe as in mandelbrot1963 writes it:
and following Black1976aa volatility tends to rise with "bad" news and fall with "good" news:
To account for those properties in our residuals, we follow Nelson:1991aa for the model notations.
Let $y_t$ denote (demeaned) asset returns. A standard decomposition writes
where $\mu_t$ captures predictable dynamics in the conditional mean (often an ARMA$(p,q)$), while $\sigma_t^2$ is the conditional variance (conditional on past/selected information: $\sigma_t^2=\sigma^2\left(\mathcal{F}_{t-1}\right)$).
The ARCH$(m)$ model Engle:1982aa specifies the conditional variance as a function of past squared innovations\footnote{See ARCH_vansteenberghe.py for the implementation.} :
The ARCH$(1)$ case is the special case
\paragraph{Unconditional variance (ARCH(1)).} Using $\mathbb{E}[z_t^2]=1$ and $\varepsilon_t=\sigma_t z_t$,
If $0 \le \alpha_1 < 1$, the stationary solution yields
which is constant over time even though $\sigma_t^2=\mathbb{V}(\varepsilon_t\mid\mathcal{F}_{t-1})$ varies with past shocks.
ARCH effects are tested on the residuals from the model. Let $\hat{\varepsilon}_t$ denote those residuals. Under the null of no ARCH, the conditional variance is constant: \[ H_0:\ \alpha_1=\cdots=\alpha_m=0 \quad\Longleftrightarrow\quad \mathbb{E}[\hat{\varepsilon}_t^2\mid\mathcal{F}_{t-1}] = \alpha_0. \] Two common diagnostics are:
\paragraph{(i) Portmanteau test on squared residuals.} If $\{\hat{\varepsilon}_t\}$ is i.i.d., then $\{\hat{\varepsilon}_t^2\}$ should be serially uncorrelated. One can therefore apply Ljung--Box tests to $\hat{\varepsilon}_t^2$ (Section (ref)).
\paragraph{(ii) Engle's ARCH LM test (recommended).} Run the auxiliary regression
and test $H_0:\alpha_1=\cdots=\alpha_m=0$.
Let $R^2$ be the coefficient of determination from (ref) and let $T$ be the sample size used in the regression. The LM statistic is
A small $p$-value indicates rejection of $H_0$ and supports the presence of ARCH effects, motivating an ARCH/GARCH-type variance specification.
Choice of the order $m$. In practice, $m$ should be large enough to capture potential short-run volatility clustering, but not so large as to lose power. Several pragmatic approaches are common:
In empirical work, it is good practice to report results for several values of $m$ and verify that conclusions about the presence of ARCH effects are not driven by a particular lag choice.
\paragraph{(iii) Li--McLeod test (robust portmanteau for heavy-tailed samples).}
In finite samples with heavy-tailed innovations, the Ljung--Box statistic applied to $\hat{\varepsilon}_t^2$ may exhibit size distortions. Li:1981aa propose a modified portmanteau statistic that remains valid under weaker moment conditions and performs better in small samples.
Let $\hat{\rho}_i$ denote the sample cross-correlation matrix at lag $i$ of the (possibly multivariate) residual process, and define \[ b_i = \mathrm{vec}(\hat{\rho}_i'). \] Let $\hat{\rho}_0$ be the sample contemporaneous correlation matrix and $\otimes$ denote the Kronecker product. The Li--McLeod statistic is
which under $H_0$ (no serial dependence up to lag $m$) satisfies \[ \mathrm{LMc} \ \overset{H_0}{\approx}\ \chi^2(k^2 m), \] where $k$ is the dimension of the residual vector.
Univariate case. When $k=1$, the statistic simplifies substantially. Let $r_i$ denote the sample autocorrelation of $\hat{\varepsilon}_t^2$ at lag $i$. Then \[ \mathrm{LMc} = T^2 \sum_{i=1}^m \frac{r_i^2}{T-i} \ \ \overset{H_0}{\approx}\ \chi^2(m). \]
This test can be viewed as a heavy-tail–robust portmanteau test on squared residuals and is particularly useful in financial applications where excess kurtosis is common.
We use the notation in Nelson:1991aa. Let $\xi_t$ denote the one-step-ahead innovation (prediction error) from the conditional mean equation, and let $\sigma_t^2$ denote its conditional variance given the information set $\mathcal{F}_{t-1}$: \[ \mathbb{E}[\xi_t \mid \mathcal{F}_{t-1}] = 0, \qquad \mathbb{V}(\xi_t \mid \mathcal{F}_{t-1}) = \sigma_t^2. \] Assume the standardized innovation is $z_t$:
where $z_t$ is often taken Gaussian, Student-$t$, or GED (Generalized Error Distribution) to allow for heavy tails.
The GARCH$(p,q)$ model generalizes ARCH by letting the conditional variance depend on both past squared innovations and past conditional variances:
with parameter restrictions \[ \omega>0,\qquad \alpha_j\ge 0,\qquad \beta_i\ge 0. \] (Equivalently, since $\xi_{t-j}^2=\sigma_{t-j}^2 z_{t-j}^2$, the ARCH term can be written as $\sum_{j=1}^{p} \alpha_j z_{t-j}^2\sigma_{t-j}^2$.)
\paragraph{Interpretation.} The $\alpha_j$ terms capture the short-run impact of past shocks on volatility (news effect), while the $\beta_i$ terms govern volatility persistence through variance feedback. The most common specification is GARCH$(1,1)$:
\paragraph{Second-order stationarity (GARCH(1,1)).} If $\alpha_1+\beta_1<1$, then $\{\xi_t\}$ has finite unconditional variance and
When $\alpha_1+\beta_1$ is close to one, volatility is highly persistent, producing long clusters of high or low variance.
\paragraph{Remark on distributions.} Assuming a non-Gaussian $z_t$ (Student-$t$, GED) typically affects likelihood-based inference and tail risk measures, but not the conditional variance recursion (ref).
For the GARCH-M model, we add the following equation $$R_t=a + b \sigma_t^2+\xi_t$$
Nelson:1991aa details how to use maximum likelihood methods to fit an ARMA(p,q) (and potentially exponential) ARCH model and discuss the regularity conditions to be checked, but as in this paper, we will not check this as this stage ans assume that the ML estimator is consistent and asymptotically normal.
We fit an AR(1)-ARCH(1) $\sigma_t^2 = \omega + \alpha_1 z_{t-1}^2\sigma_{t-1}^2$ and an AR(1)-GARCH(1,1) $\sigma_t^2 = \omega + \beta_i \sigma_{t-1}^2 + \alpha_1 z_{t-1}^2\sigma_{t-1}^2$ model on the financial index daily returns.
EGARCH models were introduced by Nelson:1991aa to accommodate three empirical features often observed in returns: (i) leverage/asymmetry (negative shocks increase future volatility more than positive shocks), (ii) fewer parameter restrictions than standard GARCH (the variance is positive by construction), and (iii) a more transparent separation between size and sign effects of shocks.
Let $R_t$ denote returns. Consider an AR(1) conditional mean and an EGARCH$(p,q)$ conditional variance:
\paragraph{Normal benchmark.} If $z_t\sim\mathcal{N}(0,1)$, then $\mathbb{E}|z_t|=\sqrt{2/\pi}$, hence the centering term $\big(|z_{t-i}|-\sqrt{2/\pi}\big)$ ensures that the intercept $\omega$ controls the unconditional level of log-variance.
The EGARCH decomposition is convenient:
Because the model is written for $\log(\sigma_t^2)$, $\sigma_t^2=\exp(\log(\sigma_t^2))>0$ automatically, so nonnegativity constraints on $\alpha_i,\beta_k$ are not required (unlike standard GARCH).
The most common specification is EGARCH$(1,1)$:
Some formulations allow $R_t$ to depend on $\sigma_t^2$ (a risk--return effect). For lecture-note simplicity and to focus on volatility dynamics, we set $c=0$ and estimate an AR(1)--EGARCH model as in (ref). Model orders $(p,q)$ can be selected using information criteria (AIC/BIC) and checked via residual diagnostics on standardized innovations $\hat{z}_t=\hat{\xi}_t/\hat{\sigma}_t$.
Nelson:1991aa suggest the EGARCH model:
as in Nelson:1991aa, we keep $c=0$ and we would need to select the order $p$, $j$, and $q$ with information criterion.
\stepcounter{question}
In equation (ref), we have a linear form of "GARCH-in-mean" with the constant $c$. We could also include any type of function (possibly non-linear) $g()$ and update equation (ref) with $R_t = a + bR_{t-1} + g\left(\sigma_t\right) + \xi_t$. We will explore section (ref) how with the package rugarch we can use the options in "mean.model" with "archm Whether to include ARCH volatility in the mean regression" and "archpow Indicates whether to use st.deviation (1) or variance (2) in the ARCH in mean regression".
Beyond checking coefficients p-values, there is the Lagrange Multiplier Test (LM Test), for a parametric model with true parameters $\theta_0$, $H_0:R\theta_0=r$, where we can test that some coefficients in $\theta_0$ are null (imposing some $0$s in $R$ and $r$).
Following Francq:2009tv, use the Wald, LM or Likelihood Ratio statistics.
Ljung-Box test on the residuals (auto-correlation), squared residuals (ARCH) and ARCH LM test (presence of ARCH).
To test the second order stationarity, you want to test $H_0: \sum_i \beta_i + \sum_j \alpha_j < 1$, with the parameters taken from equation (ref), against $H_1: \sum_i \beta_i + \sum_j \alpha_j \geq 1$.
In many financial and macroeconomic time series, the conditional mean may be relatively stable, while the conditional variance exhibits substantial fluctuations over time. Periods of calm alternate with periods of turbulence, a phenomenon commonly referred to as volatility clustering. Models with time-varying second-order moments provide a natural framework to formalize and measure such uncertainty. Uncertainty is modeled as the conditional variance of shocks. In Autoregressive Conditional Heteroskedasticity-type models, uncertainty is not an exogenous index but an endogenous, model-implied object:
where $\mathcal{F}_{t-1}$ is the information set available at time $t-1$.
Thus, ARCH/GARCH models deliver a real-time measure of perceived risk based solely on past realizations. This makes them particularly useful for:
\paragraph{Application to Inflation Uncertainty.}
While ARCH-type models were initially developed for financial returns, they are equally relevant for macroeconomic variables such as inflation. Inflation dynamics typically exhibit persistence in the mean and episodes of heightened volatility (e.g.\ oil shocks, regime shifts, pandemic episodes). ARCH-type models have been used to measure inflation uncertainty, as in grier2000effects. A natural specification therefore combines:
Let $\pi_t$ denote (annualized) inflation. A standard specification is:
where $\mathcal{D}$ may be Gaussian or Student-$t$ to allow for fat tails.
Inflation uncertainty is then defined as:
A GARCH(1,1) specification,
captures both the immediate impact of inflation surprises ($\alpha$) and the persistence of inflation uncertainty ($\beta$).
\paragraph{Economic Interpretation.}
In this framework, inflation uncertainty is the conditional variance of unexpected inflation. It reflects how difficult inflation is to forecast given available information. Importantly:
\paragraph{Mean–Variance Interactions.}
A central question in macroeconomics concerns the relationship between the level of inflation and its uncertainty. Two influential hypotheses are:
Empirically, this can be tested by augmenting the mean equation with conditional variance:
where $\gamma$ captures the effect of uncertainty on inflation dynamics.
\paragraph{Distributional Considerations.}
Inflation shocks may exhibit fat tails during crisis periods. Replacing the normal distribution with a Student-$t$ distribution often improves model fit and produces more robust measures of uncertainty:
with degrees of freedom $\nu$ estimated jointly with the volatility process.
Before\footnote{We suggest the companion script vansteenberghe_cointegration_ecm_regression.py.} applying any regression model or computing correlations, it is essential to visualize the data. In time-series analysis, plotting the levels and their transformations is the first diagnostic step: it allows us to detect trends, breaks, seasonality, outliers, and potential non-stationarity.
In this section, we refer to the second variable as “GDP” for notational simplicity. However, the series actually used in the code corresponds to:
published by \href{https://www.insee.fr/fr/statistiques/serie/001689203}{INSEE, series 001689203}. This variable is not equivalent to GDP. It measures the aggregate wage bill (in current prices, seasonally adjusted), whereas GDP captures total value added in the economy. The two series are economically related but conceptually distinct.
We use this wage-bill series for educational purposes only, in order to illustrate the mechanics of regression, correlation, unit-root testing, and cointegration in a macroeconomic context.
We perform the same data cleaning operations as for the population CSV file. Since the wage-bill series is quarterly, we align frequencies by converting the monthly population series to quarterly averages before merging the two datasets.
We then concatenate the two time series using the date index (see the \href{http://pandas.pydata.org/pandas-docs/stable/merging.html}{pandas merging documentation}) and plot them, using a secondary $y$-axis for population.
You should obtain:
We work with quarterly growth rates: \[ g_t = \frac{X_t - X_{t-1}}{X_{t-1}} \qquad \text{or equivalently} \qquad g_t \approx \Delta \log X_t . \]
We compute the quarterly changes and expect a relatively low contemporaneous correlation between wage-bill growth (labeled “GDP” in the code) and population growth.
We then inspect a scatter plot of wage-bill growth against population growth:
Two observations clearly appear as outliers: one associated with a strongly negative wage-bill growth rate (macro shock) and another with an unusually high population growth rate. We discuss formal treatment of outliers in Section (ref).
Following the methodology detailed in Section (ref), we compute the sample correlation between the two growth rates. We obtain a relatively low correlation of approximately $18\%$.
This modest correlation illustrates an important empirical lesson: even variables that are economically related in levels (e.g.\ activity and demographics) may display weak short-run comovement in growth rates.
We test "manually" the GDP for a unit root following the test described in section (ref).
We test that both the French population and GDP time series are I(1).
After we differenced the variables, we cannot check for their long term relation in level, when differenced the variables are only compared in short term models.
With two I(1) series\footnote{we define by I(n) integrated series of order n} $y_t$ and $x_t$, level comparison via linear regression would be difficult because subject to spurious regression. But two I(1) series are said to be cointegrated, if there exist a linear combination\footnote{called long-run equilibrium relationship} of them (e.g. $\beta_0 + \beta_y y_t + \beta_x x_t$) that is I(0). This means that $y_t$ and $x_t$ share a common stochastic trend and the two time series, after having applied some coefficients, cannot drift too far apart from one another in the long-run\footnote{unless the "forces" that kept them together in the past do not hold any longer}. If there exist\footnote{$\beta_y$ and $\beta_z$ are not unique (one can add a constant and the relationship would still hold), also note that $\beta_0$ can be null and the following still applies} $\beta_y$ and $\beta_x$ such that $E(\beta_0 + \beta_y y_t + \beta_x x_t) = 0$, then there is a long term relationship between $y_t$ and $x_t$ and one can test the following model:
and from our statements above, one can deduct that the cointegrating residuals $\epsilon_t$ have to be $I(0)$ if $y_t$ and $x_t$ are cointegrated\footnote{but as demonstrated in Phillips and Ouliaris (1990), one cannot use the traditional Augmented Dickey-Fuller for this test because of the spurious nature of the regression, furthermore, the the t-statistics for $\beta$ will not be asymptotically normal if $\epsilon_t$ are serially correlated}.
We can then test for cointegration between variable with a cointegration test: the null hypothesis is no cointegration. The second output give in the p-value of the test.
For our variables, we do not reject the null hypothesis and assume no cointegrations.
We suggest to build two cointegrated time series with the system of equations:
with both $\epsilon_t$ and $\nu_t$ white noises.
We can show that both series are I(1) and cointegrated. Visually the link might not be obvious:
Once we applied the cointegrating regression and found the coefficients, we can show the "forces de rappel" where there seem to be some forces that brings $y_t$ to "follow" $x_t$:
We can then study the following Error Correction Model, with a grain of salt\footnote{e.g. "ECM system may be estimated by seemingly unrelated regressions (SUR) to increase efficiency if the number of lags in the two equations are different"} with regard to the consistency of the estimation by OLS as detailed in zivot2003:
We take the Euribor 1 year (FM.M.U2.EUR.RT.MM.EURIBOR1YD_.HSTA) and the Euribor 3 months (FM.M.U2.EUR.RT.MM.EURIBOR3MD_.HSTA) from the ECB statistical warehouse, check for cointegration between the two time series and build an ECM if relevant.
Noting that both time series are I(1) and cointegrated, there is no unit root in the error terms $\epsilon_t$ of the regression: $$Y_t = \alpha + \beta X_t + \epsilon_t$$
Remember that if our model is correct, then $E(\epsilon_t) = 0$
We first estimate the long run relationship from the spurious model: $$Y_t = \alpha + \beta X_t + \epsilon_t$$
and we estimate the disequilibrium errors or the cointegrating residuals: $$\hat{\epsilon}_t = Y_t - \hat{\alpha} - \hat{\beta} X_t$$
We then compute the error correction model were all variable are now stationary: $$\Delta Y_t = \gamma_0 + \gamma_1 \Delta X_t + \gamma_2 \hat{\epsilon}_{t-1} + \nu_t$$
We can say that our model explains $78\%$ of the variance of the observed changes. And that for a $1\%$ increase in the Eur