# LARCH: Large Language Model-based Automatic Readme Creation with Heuristics

Yuta Koreeda<EMAIL_ADDRESS>Hitachi, Ltd., Research and
Development Group1-280 HigashikoigakuboKokubunjiTokyoJapan185-0014 ,  Terufumi
Morishita<EMAIL_ADDRESS>Hitachi, Ltd., Research and
Development Group1-280 HigashikoigakuboKokubunjiTokyoJapan185-0014 ,  Osamu
Imaichi<EMAIL_ADDRESS>Hitachi, Ltd., Research and Development
Group1-280 HigashikoigakuboKokubunjiTokyoJapan185-0014 and  Yasuhiro Sogawa
<EMAIL_ADDRESS>Hitachi, Ltd., Research and Development
Group1-280 HigashikoigakuboKokubunjiTokyoJapan185-0014

(2023)

###### Abstract.

Writing a _readme_ is a crucial aspect of software development as it plays a
vital role in managing and reusing program code. Though it is a pain point for
many developers, automatically creating one remains a challenge even with the
recent advancements in large language models (LLMs), because it requires
generating an abstract description from thousands of lines of code. In this
demo paper, we show that LLMs are capable of generating a coherent and
factually correct readmes if we can identify a code fragment that is
representative of the repository. Building upon this finding, we developed
LARCH (LLM-based Automatic Readme Creation with Heuristics) which leverages
representative code identification with heuristics and weak supervision.
Through human and automated evaluations, we illustrate that LARCH can generate
coherent and factually correct readmes in the majority of cases, outperforming
a baseline that does not rely on representative code identification. We have
made LARCH open-source and provided a cross-platform Visual Studio Code
interface and command-line interface, accessible at
https://github.com/hitachi-nlp/larch. A demo video showcasing LARCH’s
capabilities is available at https://youtu.be/ZUKkh5ED-O4.

large language model, software development, weak supervision

††journalyear: 2023††copyright: acmlicensed††conference: Proceedings of the
32nd ACM International Conference on Information and Knowledge Management;
October 21–25, 2023; Birmingham, United Kingdom††booktitle: Proceedings of the
32nd ACM International Conference on Information and Knowledge Management
(CIKM ’23), October 21–25, 2023, Birmingham, United Kingdom††price:
15.00††doi: 10.1145/3583780.3614744††isbn: 979-8-4007-0124-5/23/10††ccs:
Software and its engineering Documentation††ccs: Computing methodologies
Natural language generation

Figure 1. The overview of LARCH (Large language model-based Automatic Readme
Creation with Heuristics)

The overview of LARCH (large language model-based automatic readme creation
with heuristics). The figure starts from a screenshot of Visual Studio Code in
which the plugin is being brought up. The whole code repository consisting of
hundreds of files are being transferred to the representative code
identification module. Then, a single file of representative code is extracted
from the repository and sent to the prompt creation module. Finally, the
prompt is sent to a large language model and the outcome is displayed on
Visual Studio Code.

##  1\. Introduction

Recent advances in AI, especially large language models (LLMs) (Bommasani et
al., 2021), are revolutionalizing software development through code search
(Husain et al., 2019), program repair (Yasunaga and Liang, 2020), code
generation (Chen et al., 2021) and many other applications. However, assisting
developers with documentations, which are as important as code itself, is not
adequately addressed even though it is a pain point for many developers (de
Souza et al., 2005). In particular, assisting developers write _readme_ s is
utmost important as it is the most written form of documentation and having no
readme essentially makes code unreusable.

The state-of-the-art in assisting developers write a readme is by merely
presenting them with a static template or populating the template based on
user input, but these solutions do not actually help them write its content.
Previous works have shown that LLMs can generate class-/function-level code
comments (Wang et al., 2021; Chen et al., 2021). However, comment generation
merely involves generating a concrete comment from dozens of lines of code.
Generating a readme remains a challenge as it requires generating an abstract
summarization of thousands or even millions of lines of code.

In this demo paper, we show that LLMs are capable of generating a coherent and
factually correct readme if we can identify a code fragment that is
_representative_ (gives overview) of the repository. Based on this finding, we
developed LARCH (LLM-based Automatic Readme Creation with Heuristics) which is
based on representative code identification with heuristics and weak
supervision (Figure 1). Our contributions are as follows:

  * •

We developed LARCH, the first system to generate coherent and factually
correct readmes utilizing LLMs.

  * •

We show the efficacy of our approach through both human and automated
evaluation.

  * •

We implemented and open sourced111https://github.com/hitachi-nlp/larch LARCH
along with a cross-platform Visual Studio Code (VSCode) interface and command
line interface (CLI).

In our demo, attendees will have chance to test our system against code of
their choice. You can find the demo video at https://youtu.be/ZUKkh5ED-O4.

##  2\. Retrieval-Augmented Language Models

A language model is a probability distribution over sequences of tokens and it
can generate a token sequence by iteratively calculating probabilities of
$(i+1)$-th token given the _context_ of preceding $i$ tokens
$\\{x_{j}\\}_{j\leq i}$ that were already generated. In order to generate
coherent readme for each repository, we need to carry out generation
contextualized by the repository information $\mathcal{R}$. Following the
recent _prompting_ paradigm (Liu et al., 2023b), we feed the repository
information $\mathcal{R}$ as sequences of tokens $\\{r_{i}\\}$:

(1)
$p(x_{0},\dotsc,x_{n}|\mathcal{R})=\prod\nolimits_{i=0}^{n-1}p(x_{i+1}|r_{0},\dotsc,r_{|\mathcal{R}|},x_{0},\dotsc,x_{i}).$

Recently dominant Transformer-based LMs (Vaswani et al., 2017) suffer from a
quadric computational complexity against the the sequence lengths. Since
viable context lengths of existing models are much shorter than the average
repository size, we need to summarize an input repository to a fixed-size
token sequence. We follow _retrieval-augmented language models_ approach (Shi
et al., 2023) and retrieve necessary information from each repository and
insert its tokens to $\mathcal{R}$. Through pilot studies, we found that
identifying representative code is the key to readme generation (as
demonstrated in Section 4), hence we developed an representative code
identification method as described in Section 3.1.

Table 1. The labeling functions (LFs) and the features for the entry point identification. | Values for
---|---
Description | LF† | Feature
_File content_
| 1. | Contains a string “main” in a function name | 1, 0 | 1, 0
| 2. | Contains an argument parser | 1, 0 | 1, 0
| 3. | Contains a web framework (such as Flask) | 1, 0 | 1, 0
| 4a. | Too short ($<200$ characters) | -1, 0 | —
| 4b. | Content length (# characters) | — | int
_Directory information_
| 5. | Contains substring “main” in the file name | 1, 0 | 1, 0
| 6. | Has entry point-ish name (such as “cli.py”) | 1, 0 | 1, 0
| 7. | Is “__init__.py” | -1, 0 | 1, 0
| 8. | Has a test-ish name (i.e., starts with “test_”) | -1, 0 | 1, 0
| 9. | Directory depth from the project root | — | int
_Static code analysis_
| 10a. | Is the top of the import tree | 1, 0 | —
| 10b. | Distance from the top in the import tree | — | int
| 11. | Is the bottom of the import tree | -1, 0 | 1, 0
| 12. | # imports (within the repository) | — | int
| 13. | # importers (within the repository) | — | int
| 14a. | Contains a class inherited $\geq 3$ times | 1, 0 | —
| 14b. | # classes inheriting a class from this file | — | int
_Oracle_
| 15. | Has the same file name as the repository | 1, 0 | —
| 16. | Listed as “entry point” in “setup.py” | 1, -1, 0‡ | —
| 17. | Imported in the reference readme | 1, -1, 0‡ | —
†1 if it is likely to be a representative code file, -1 if not, and 0 if it
abstains. ‡0 if it is not 1 and there exists at least one file in the same
repository that is 1.

##  3\. LARCH: LLM-based Automatic Readme Creation with Heuristics

The overview of LARCH is shown in Figure 1. Users can launch LARCH from
VSCode. LARCH aggregates code repository from the current workspace, which is
sent to the API server. LARCH identifies the most representative code of the
repository using heuristics-based features and gradient boosting trees
(Section 3.1). Then, a prompt is constructed from the extracted code (Section
3.2) and it is sent to either an external LLM API or a local LLM (Section
3.3). Finally, the generated readme is sent back to the user and shown on the
editor. The whole process takes about 20 seconds.

The API server comes with CLI which can communicate with API server or
generate a readme without launching the server. They are distributed as a
Python pip package and hence cross-platform. The VSCode interface is
distributed as a VSCode plugin (a “.vsix” file) and is also cross-platform.

In this paper, we focus on Python projects as it is the most prominent
programming language used in the machine learning community. Nevertheless, our
framework can be extended to other languages as well.

###  3.1. Weak Supervision for Representative Code Identification

Representative code can take many forms; It can be an entry point of an
application, a facade in facade design pattern (Gamma et al., 1994), or a base
class in object-oriented libraries. Since such complex concept cannot be
captured by static program analysis, we propose to employ heuristics-based
features that consider diverse properties of a repository, and use machine
learning to identify representative code file222It can be a different
granularity such as functions, but we chose files as it is simple and their
lengths match LLMs’ context lengths..

The representative code identification problem has properties that  (1)
annotation is quite costly as it requires careful inspection of each
repository, and  (2) we can obtain lots of unlabeled public repositories  .
Hence, we decided to take data programming paradigm (Ratner et al., 2017), a
weak supervision approach, where we handcraft a set of heuristics to create
silver labels for training a machine learning model. More specifically, we
implemented labeling functions where $j$-th function takes $i$-th file and
returns a noisy label ($\in\\{-1,1\\}$) or abstains ($=0$;
$\Lambda_{i,j}\in\\{-1,1,0\\}$; $\mathbf{\Lambda}=\\{\Lambda_{i,j}\\}$). We
can then recover the accuracies of these labeling functions and label
posterior for $i$-th file $p(y_{i}|\mathbf{\Lambda})$ (where
$y_{i}\in\\{-1,1\\}$) by solving a matrix completion-style problem333We used
Snorkel (https://www.snorkel.org/) (Ratner et al., 2019)..

We developed 14 labeling functions (Table 1). Notably, we utilize oracle
information from reference readmes. This wouldn’t be available in an ordinary
weak supervision setting but it should serve as a strong cue for identifying
representative code.

Figure 2. The prompt design

The figure shows the prompt template. It starts with a phrase “Here is the
entrypoint of a Python project “larch”:” which is followed by the extracted
representative code. It is followed by text “This program has following
files:” and ten randomly selected file names. Finally, it contains a phrase
“Write a detailed readme in markdown:” to mark the start of generation.

After obtaining $\\{p(y_{i}=1|\mathbf{\Lambda})\\}$, we train gradient
boosting trees (Chen and Guestrin, 2016) to identify representative code
files. We use 14 features that are similar to our labeling functions (Table
1). We omit oracle information that wouldn’t be available at inference time,
and we stop discretizing features to Boolean as we have more flexibility in
the values. Instead of formulating the problem as file-wise binary
classification, we formulated it as a learning-to-rank problem of files within
each repository. This formulation is more appropriate as our objective is to
pick a single file from each repository. This trick can also improve the
overall accuracy by removing repository-level biases of labeling functions
(e.g., a repository may contain many files whose names contain “main”).

###  3.2. Prompt Design

Previous studies have shown that the design of prompt has a significant effect
on downstream tasks (Reynolds and McDonell, 2021; Kojima et al., 2022). Hence
we carefully designed a prompt template that utilizes extracted representative
code (Figure 2). We found that, at least for GPT-3 (Ouyang et al., 2022),
imperative sentences are better than declarative sentences (e.g., “I wrote a
readme for this program:”). Specifying “markdown” and “Python” did not help
much in most cases, but it seems to avoid catastrophic mistakes especially
when the input code is short. If the project name is not given, LLMs generally
come up with a name that most matches the code, so we made project name an
optional input. Having file names helps for projects that implement many
variations of a single functionality.

###  3.3. Large Language Models

We utilized OpenAI API’s GPT-3 “davinci-text-003”444https://openai.com/api/
for the LLM. We utilized prompt length of 3,000 tokens and maximum generation
length of 910 tokens.

Table 2. Evaluation of readme generation with representative code
identification and the random file baseline

(a) Overall human evaluation Context |  Useless |  Fair |  Useful
---|---|---|---
Random file |  40% |  20% |  40%
Representative code (ours) |  15% |  20% |  65%
“Useful” when it can be adopted with small fixes, “Fair” when it may be useful
as a reference, and “Useless” otherwise.

(b) Fine-grained human evaluation (% of positive assessments) | Context
---|---
Criteria | Random file | Representative code (ours)
Includes project goal | 100% | 100%
Includes instruction | 100% | 100%
Grammatical correctness† | 100% | 100%
Markdown correctness† | 100% | 100%
Factual correctness (text)‡ | 55% | 75%
Factual correctness (code)‡ | 30% | 65%
†Percentage of “Good” from the choices of “Bad”, “Fair” and “Good”. ‡“Good”
when it is mostly correct (e.g., we don’t require an example code to run as
is, as in a human-written readme), “Bad” otherwise.

(c) Automatic evaluation | ROUGE score $\uparrow$
---|---
Context |  1 |  2 |  L
Random file |  20.9 |  4.9 |  10.7
Representative code (ours) |  22.0 |  5.9 |  11.4
Each score considers different $n$-grams. See (Lin, 2004) for details.

##  4\. Experiments

We evaluated LARCH with both human and automatic evaluation. Since previous
studies are limited, we compare LARCH against a baseline that uses a randomly
selected Python file instead of representative code, as we believe that
representative code identification is the key component of LARCH.

(a) vitsalis/PyCG

(b) raphaelquast/EOmaps

Figure 3. Excerpts of the generated readmes

These set of figures show samples of generated readmes. Both examples show
that generated readmes are coherent and mostly factually correct. However, it
also points out that LARCH made a minor mistake of writing GeoTIFF as GeoJSON.
It also shows that it gets the license term completely incorrect in one
example; It says that the repository is MIT but it is actually GPL.

We collected public repositories from GitHub555https://github.com/, each of
which  (1) has more than 100 stars,  (2) written in Python,  (3) is smaller
than 500 MB in size with less than 1,000 files,  (4) contains English readme
in markdown format, and  (5) created after the publication of GPT-3 (Brown et
al., 2020) (11th June 2020)  . We randomly sampled 1,500 repositories for
automatic evaluation, 20 of which were also used for human evaluations. We
removed a readme and “setup.py” from each repository and kept readme as
reference data. Each repository contained, on average, 38.3 files and 103,302
tokens.

For the human evaluation, we compared the two systems in a double-blind
setting (pairs of outputs were anonymized and presented in a random order).
The results of the human evaluation are shown in Table 2. For all the
repositories, both systems managed cover both project goals and instructions
(“getting started”) and they were coherent both in terms of grammar and
markdown formatting. LARCH, however, performed much better in terms of factual
correctness. This lead to the significant improvement in “overall usefulness”
($p=0.009$, Wilcoxon signed-rank test) — our system performed equally or
better than the baseline in 95% of the repositories.

For the automatic evaluation, we compared ROUGE score (Lin, 2004) of generated
readmes to the reference readmes, a common metric used in summarization. As
shown in Table 2(c), LARCH outperformed the baseline in all metrics.

A couple of example outputs are shown in Figure 3. Figure 3(a) is a
straightforward application and LARCH correctly identified the project goal
and its example usage. LARCH chose a file with the entry point where an
argument parser and main operations are located. Figure 3(b) is a class-based
library and LARCH chose a file with the main class container. LARCH got the
project features and its usage mostly correct, but got the license completely
wrong. It is natural as we did not incorporate license information to the
prompt. This is an easily amendable problem with some engineerings and we
leave it for the future work.

##  5\. Related Works

As discussed in Section 1, prevailing approaches to aiding developers in
crafting informative readmes predominantly rely on template-based methods.
However, since and around the submission of this paper on June 16, 2023,
multiple relevant works have emerged that warrant a discussion.

StarCoder (Li et al., 2023) (released on May 4) is the 1.5B parameters state-
of-the-art model for code generation. While readme generation is more of
natural language generation than code generation, we believe it is worth
running experiments with StarCoder and other code generation models (Chen et
al., 2021) in the future work. There are also multiple new models that support
much longer context lengths; GPT-4-32K (OpenAI, 2023b, a) (released on July 6)
supports 32K tokens and Claude-2 (Anthropic, 2023) (released on July 11)
supports 100K. These may lessen the needs for our representative code
identification, but we argue that it remains important because  (1) there
exist many repositories that still do not fit onto these larger context
lengths,  (2) the performance of LLMs tend to degrade for lengthy inputs even
if their positional encodings _support_ them (Liu et al., 2023a), and  (3)
processing long inputs requires more compute  . We would nevertheless like to
compare and incorporate newer LLMs to our work in the future.

README-AI666https://github.com/eli64s/readme-ai/tree/v0.0.5 (version 0.0.1
released on June 28) is a Python library that generates readmes with a LLM but
with a slightly different approach. README-AI summarizes each source file into
a short description first, and generates a readme from the concatenated
descriptions. This makes each source file significantly shorter and allows
fitting more than one source file to the context. Nevertheless, we believe
that the same arguments as the long context LLMs apply to READM-AI; it does
not completely solve the context length limit and it can actually be used
along side with our method. In the future work, we would like to
quantitatively compare two approaches and evaluate how they perform if they
are put together.

##  6\. Conclusion

As presented in Section 4, LARCH can generate coherent and factually correct
readme in majority of cases. We have shown that our representative code
identification approach yields much better generation than the baseline. While
there exist risks that LARCH may generate factually incorrect readme,
developers can always fix the result. Since LARCH is straightforward to use
and reading readme is much easier than writing one, LARCH can assist
developers without having negative effect to the community.

For future work, we will extend our framework to different programming
languages.

## Ethical Consideration

While LARCH significantly improves factual correctness from the baseline, it
can still get facts wrong as demonstrated in Section 4. Misinformation in
readmes can have negative effect to the users as it may introduce bugs or
causes legal issues with regard to the licencing. That being said, we intend
LARCH to be used by developers of a repository themselves, hence they can
always neglect or fix the result. Since LARCH is straightforward to use and
reading readme is much easier than writing one, LARCH can still assist
developers without having negative effect to the community.

While we did not run computationally expensive pretraining of LLMs, LARCH
still relies on a LLM at the inference which has an unignorable carbon
footprint. Yet, we believe this energy consumption can be justified by
resources that we can potentially save by assisting developers.

###### Acknowledgements.

We used computational resource of AI Bridging Cloud Infrastructure (ABCI)
provided by the National Institute of Advanced Industrial Science and
Technology (AIST) for the experiments. We would like to thank Dr. Masaaki
Shimizu for arranging the computational environment and Takuo Shigetani for
helping us with the UI implementation.

## References

  * (1)
  * Anthropic (2023) Anthropic. 2023\.  Model Card and Evaluations for Claude Models.  (2023).  https://www-files.anthropic.com/production/images/Model-Card-Claude-2.pdf Accessed August 16, 2023.
  * Bommasani et al. (2021) Rishi Bommasani, Drew A. Hudson, Ehsan Adeli, Russ Altman, Simran Arora, Sydney von Arx, Michael S. Bernstein, Jeannette Bohg, Antoine Bosselut, Emma Brunskill, Erik Brynjolfsson, S. Buch, Dallas Card, Rodrigo Castellon, Niladri S. Chatterji, Annie S. Chen, Kathleen A. Creel, Jared Davis, Dora Demszky, Chris Donahue, Moussa Doumbouya, Esin Durmus, Stefano Ermon, John Etchemendy, Kawin Ethayarajh, Li Fei-Fei, Chelsea Finn, Trevor Gale, Lauren E. Gillespie, Karan Goel, Noah D. Goodman, Shelby Grossman, Neel Guha, Tatsunori Hashimoto, Peter Henderson, John Hewitt, Daniel E. Ho, Jenny Hong, Kyle Hsu, Jing Huang, Thomas F. Icard, Saahil Jain, Dan Jurafsky, Pratyusha Kalluri, Siddharth Karamcheti, Geoff Keeling, Fereshte Khani, O. Khattab, Pang Wei Koh, Mark S. Krass, Ranjay Krishna, Rohith Kuditipudi, Ananya Kumar, Faisal Ladhak, Mina Lee, Tony Lee, Jure Leskovec, Isabelle Levent, Xiang Lisa Li, Xuechen Li, Tengyu Ma, Ali Malik, Christopher D. Manning, Suvir P. Mirchandani, Eric Mitchell, Zanele Munyikwa, Suraj Nair, Avanika Narayan, Deepak Narayanan, Benjamin Newman, Allen Nie, Juan Carlos Niebles, Hamed Nilforoshan, J. F. Nyarko, Giray Ogut, Laurel Orr, Isabel Papadimitriou, Joon Sung Park, Chris Piech, Eva Portelance, Christopher Potts, Aditi Raghunathan, Robert Reich, Hongyu Ren, Frieda Rong, Yusuf H. Roohani, Camilo Ruiz, Jack Ryan, Christopher Ré, Dorsa Sadigh, Shiori Sagawa, Keshav Santhanam, Andy Shih, Krishna Parasuram Srinivasan, Alex Tamkin, Rohan Taori, Armin W. Thomas, Florian Tramèr, Rose E. Wang, William Wang, Bohan Wu, Jiajun Wu, Yuhuai Wu, Sang Michael Xie, Michihiro Yasunaga, Jiaxuan You, Matei A. Zaharia, Michael Zhang, Tianyi Zhang, Xikun Zhang, Yuhui Zhang, Lucia Zheng, Kaitlyn Zhou, and Percy Liang. 2021.  On the Opportunities and Risks of Foundation Models.  _arXiv_ (2021).  https://arxiv.org/abs/2108.07258
  * Brown et al. (2020) Tom Brown, Benjamin Mann, Nick Ryder, Melanie Subbiah, Jared D Kaplan, Prafulla Dhariwal, Arvind Neelakantan, Pranav Shyam, Girish Sastry, Amanda Askell, Sandhini Agarwal, Ariel Herbert-Voss, Gretchen Krueger, Tom Henighan, Rewon Child, Aditya Ramesh, Daniel Ziegler, Jeffrey Wu, Clemens Winter, Chris Hesse, Mark Chen, Eric Sigler, Mateusz Litwin, Scott Gray, Benjamin Chess, Jack Clark, Christopher Berner, Sam McCandlish, Alec Radford, Ilya Sutskever, and Dario Amodei. 2020.  Language Models are Few-Shot Learners. In _Advances in Neural Information Processing Systems_.  https://proceedings.neurips.cc/paper/2020/file/1457c0d6bfcb4967418bfb8ac142f64a-Paper.pdf
  * Chen et al. (2021) Mark Chen, Jerry Tworek, Heewoo Jun, Qiming Yuan, Henrique Ponde de Oliveira Pinto, Jared Kaplan, Harri Edwards, Yuri Burda, Nicholas Joseph, Greg Brockman, Alex Ray, Raul Puri, Gretchen Krueger, Michael Petrov, Heidy Khlaaf, Girish Sastry, Pamela Mishkin, Brooke Chan, Scott Gray, Nick Ryder, Mikhail Pavlov, Alethea Power, Lukasz Kaiser, Mohammad Bavarian, Clemens Winter, Philippe Tillet, Felipe Petroski Such, Dave Cummings, Matthias Plappert, Fotios Chantzis, Elizabeth Barnes, Ariel Herbert-Voss, William Hebgen Guss, Alex Nichol, Alex Paino, Nikolas Tezak, Jie Tang, Igor Babuschkin, Suchir Balaji, Shantanu Jain, William Saunders, Christopher Hesse, Andrew N. Carr, Jan Leike, Josh Achiam, Vedant Misra, Evan Morikawa, Alec Radford, Matthew Knight, Miles Brundage, Mira Murati, Katie Mayer, Peter Welinder, Bob McGrew, Dario Amodei, Sam McCandlish, Ilya Sutskever, and Wojciech Zaremba. 2021\.  Evaluating Large Language Models Trained on Code.  _arXiv_ (2021).  https://arxiv.org/abs/2107.03374
  * Chen and Guestrin (2016) Tianqi Chen and Carlos Guestrin. 2016.  XGBoost: A Scalable Tree Boosting System. In _Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining_.  https://doi.org/10.1145/2939672.2939785
  * de Souza et al. (2005) Sergio Cozzetti B. de Souza, Nicolas Anquetil, and Káthia M. de Oliveira. 2005.  A Study of the Documentation Essential to Software Maintenance. In _Proceedings of the 23rd Annual International Conference on Design of Communication_.  https://doi.org/10.1145/1085313.1085331
  * Gamma et al. (1994) Erich Gamma, Richard Helm, Ralph Johnson, and John Vlissides. 1994.  _Design Patterns: Elements of Reusable Object-Oriented Software_.  Addison-Wesley Professional.
  * Husain et al. (2019) Hamel Husain, Ho-Hsiang Wu, Tiferet Gazit, Miltiadis Allamanis, and Marc Brockschmidt. 2019.  CodeSearchNet Challenge: Evaluating the State of Semantic Code Search.  (2019).  https://arxiv.org/abs/1909.09436
  * Kojima et al. (2022) Takeshi Kojima, Shixiang Shane Gu, Machel Reid, Yutaka Matsuo, and Yusuke Iwasawa. 2022\.  Large Language Models are Zero-Shot Reasoners. In _Advances in Neural Information Processing Systems_.  https://openreview.net/forum?id=e2TBb5y0yFf
  * Li et al. (2023) Raymond Li, Loubna Ben Allal, Yangtian Zi, Niklas Muennighoff, Denis Kocetkov, Chenghao Mou, Marc Marone, Christopher Akiki, Jia Li, Jenny Chim, Qian Liu, Evgenii Zheltonozhskii, Terry Yue Zhuo, Thomas Wang, Olivier Dehaene, Mishig Davaadorj, Joel Lamy-Poirier, João Monteiro, Oleh Shliazhko, Nicolas Gontier, Nicholas Meade, Armel Zebaze, Ming-Ho Yee, Logesh Kumar Umapathi, Jian Zhu, Benjamin Lipkin, Muhtasham Oblokulov, Zhiruo Wang, Rudra Murthy, Jason Stillerman, Siva Sankalp Patel, Dmitry Abulkhanov, Marco Zocca, Manan Dey, Zhihan Zhang, Nour Fahmy, Urvashi Bhattacharyya, Wenhao Yu, Swayam Singh, Sasha Luccioni, Paulo Villegas, Maxim Kunakov, Fedor Zhdanov, Manuel Romero, Tony Lee, Nadav Timor, Jennifer Ding, Claire Schlesinger, Hailey Schoelkopf, Jan Ebert, Tri Dao, Mayank Mishra, Alex Gu, Jennifer Robinson, Carolyn Jane Anderson, Brendan Dolan-Gavitt, Danish Contractor, Siva Reddy, Daniel Fried, Dzmitry Bahdanau, Yacine Jernite, Carlos Muñoz Ferrandis, Sean Hughes, Thomas Wolf, Arjun Guha, Leandro von Werra, and Harm de Vries. 2023\.  StarCoder: May the Source be with You!  _arXiv_ (2023).  https://arxiv.org/abs/2305.06161v1
  * Lin (2004) Chin-Yew Lin. 2004\.  ROUGE: A Package for Automatic Evaluation of Summaries. In _Proceedings of the Workshop on Text Summarization Branches Out_. 74–81.  https://aclanthology.org/W04-1013
  * Liu et al. (2023a) Nelson F. Liu, Kevin Lin, John Hewitt, Ashwin Paranjape, Michele Bevilacqua, Fabio Petroni, and Percy Liang. 2023a.  Lost in the Middle: How Language Models Use Long Contexts.  _arXiv_ (2023).  https://arxiv.org/abs/2307.03172v2
  * Liu et al. (2023b) Pengfei Liu, Weizhe Yuan, Jinlan Fu, Zhengbao Jiang, Hiroaki Hayashi, and Graham Neubig. 2023b.  Pre-Train, Prompt, and Predict: A Systematic Survey of Prompting Methods in Natural Language Processing.  _Comput. Surveys_ 55, 9 (2023).  https://doi.org/10.1145/3560815
  * OpenAI (2023a) OpenAI. 2023a.  GPT-4.  _OpenAI Blog_ (2023).  https://openai.com/research/gpt-4 Accessed August 16, 2023.
  * OpenAI (2023b) OpenAI. 2023b.  GPT-4 Technical Report.  _arXiv_ (2023).  https://arxiv.org/abs/2303.08774v3
  * Ouyang et al. (2022) Long Ouyang, Jeffrey Wu, Xu Jiang, Diogo Almeida, Carroll Wainwright, Pamela Mishkin, Chong Zhang, Sandhini Agarwal, Katarina Slama, Alex Gray, John Schulman, Jacob Hilton, Fraser Kelton, Luke Miller, Maddie Simens, Amanda Askell, Peter Welinder, Paul Christiano, Jan Leike, and Ryan Lowe. 2022\.  Training Language Models to Follow Instructions with Human Feedback. In _Advances in Neural Information Processing Systems_.  https://openreview.net/forum?id=TG8KACxEON
  * Ratner et al. (2017) Alexander Ratner, Stephen H. Bach, Henry Ehrenberg, Jason Fries, Sen Wu, and Christopher Ré. 2017.  Snorkel: Rapid Training Data Creation with Weak Supervision. In _Proceedings of the VLDB Endowment_.  https://doi.org/10.14778/3157794.3157797
  * Ratner et al. (2019) Alexander Ratner, Braden Hancock, Jared Dunnmon, Frederic Sala, Shreyash Pandey, and Christopher Ré. 2019.  Training Complex Models with Multi-Task Weak Supervision. In _Proceedings of the Thirty-Third AAAI Conference on Artificial Intelligence and Thirty-First Innovative Applications of Artificial Intelligence Conference and Ninth AAAI Symposium on Educational Advances in Artificial Intelligence_.  https://doi.org/10.1609/aaai.v33i01.33014763
  * Reynolds and McDonell (2021) Laria Reynolds and Kyle McDonell. 2021.  Prompt Programming for Large Language Models: Beyond the Few-Shot Paradigm. In _Extended Abstracts of the 2021 CHI Conference on Human Factors in Computing Systems_.  https://doi.org/10.1145/3411763.3451760
  * Shi et al. (2023) Weijia Shi, Sewon Min, Michihiro Yasunaga, Minjoon Seo, Rich James, Mike Lewis, Luke Zettlemoyer, and Wen-tau Yih. 2023\.  REPLUG: Retrieval-Augmented Black-Box Language Models.  _arXiv_ (2023).  https://doi.org/10.48550/ARXIV.2301.12652
  * Vaswani et al. (2017) Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan N Gomez, Łukasz Kaiser, and Illia Polosukhin. 2017.  Attention is All you Need. In _Advances in Neural Information Processing Systems_.  https://proceedings.neurips.cc/paper/2017/file/3f5ee243547dee91fbd053c1c4a845aa-Paper.pdf
  * Wang et al. (2021) Yue Wang, Weishi Wang, Shafiq Joty, and Steven C.H. Hoi. 2021\.  CodeT5: Identifier-aware Unified Pre-trained Encoder-Decoder Models for Code Understanding and Generation. In _Proceedings of the 2021 Conference on Empirical Methods in Natural Language Processing_.  https://doi.org/10.18653/v1/2021.emnlp-main.685
  * Yasunaga and Liang (2020) Michihiro Yasunaga and Percy Liang. 2020.  Graph-Based, Self-Supervised Program Repair from Diagnostic Feedback. In _Proceedings of the 37th International Conference on Machine Learning_.  https://dl.acm.org/doi/abs/10.5555/3524938.3525939


# Application of Kalman Filter in Stochastic Differential Equations

Wencheng Bao, Shi Feng
Department of Applied Physics and Applied Mathematics, Graduate candidate
Columbia University, New York, NY
<EMAIL_ADDRESS><EMAIL_ADDRESS>&Kaiwen Zhang
Department of Applied Physics and Applied Mathematics, Undergraduate candidate
Columbia University, New York, NY
<EMAIL_ADDRESS>
###### Abstract

In areas such as finance, engineering, and science, we often face situations
that change quickly and unpredictably. These situations are tough to handle
and require special tools and methods capable of understanding and predicting
what might happen next. Stochastic Differential Equations (SDEs) are renowned
for modeling and analyzing real-world dynamical systems. However, obtaining
the parameters, boundary conditions, and closed-form solutions of SDEs can
often be challenging. In this paper, we will discuss the application of Kalman
filtering theory to SDEs, including Extended Kalman filtering and Particle
Extended Kalman filtering. We will explore how to fit existing SDE systems
through filtering and track the original SDEs by fitting the obtained closed-
form solutions. This approach aims to gather more information about these
SDEs, which could be used in various ways, such as incorporating them into
parameters of data-based SDE models.

_K_ eywords Kalman Filter $\cdot$ Extended Kalman Filter $\cdot$ Particle
Kalman Filter $\cdot$ Heston $\cdot$ Stochastic Differential Equations

##  1 Introduction

Stochastic Differential Equations (SDEs) are formed by adding one white noise
term to an ordinary differential equation. These equations are particularly
useful in modeling systems or processes that exhibit random behaviors. In real
life, SDEs have many applications, such as the Black-Scholes Model[1] used in
option pricing, the Fokker-Planck Equation[2] in physics describing the time
evolution of probability densities under random forces, and the Lotka-Volterra
Equations with Noise[3] in describing classic predator-prey relations.

However, not all SDEs have closed-form solutions. Additionally, there arises
the question of how to apply known SDEs to actual data, including how to
select the most appropriate parameters and set the appropriate boundary
conditions. These questions have been studied up to the present day. For
example, the Heston model is known as a model for option pricing [4]. To solve
parameter estimation in the Heston model, various methods have been tried,
including Maximum Likelihood Estimation (MLE) based on historical equity
prices, the Efficient Method of Moments (EMM) algorithm, the Monte Carlo
Markov Chain (MCMC) algorithm, and Quasi-Maximum Likelihood Estimation
(QMLE)[5]. Nevertheless, when it comes to capturing elusive factors like
unknown volatility, filtering theory has emerged as a pivotal approach,
offering new avenues for understanding and managing the complexities inherent
in SDEs.

The filtering theory studies the complex challenges of estimating systems, as
the estimation system itself is incomplete and susceptible to noise
interference[6]. This field is mainly applied in two key situations: first,
direct observation of variables is impractical and requires indirect
measurement; second, it is necessary to integrate multiple variables, each of
which is prone to random fluctuations. This article aims to focus on the
Kalman filtering and its various extensions, exploring their roles and
significance in the above SDEs. By applying these advanced filtering methods,
we can approximate the behavior of SDEs even in the absence of neat, closed-
form solutions, offering an approach to understanding and navigating the
stochastic dynamics of real-world phenomena.

##  2 Technique Simulation Review

On the problem of learning and fitting time-series data with uncertainties,
there is abundant literature using advanced computer tools and programs to
perform the task. These tools include various types of deep neural networks,
including multi-layer perceptron (MLP)[7], recurrent neural networks such as
gated recurrent unit (GRU)[8] and long short term memory (LSTM)[9]), and
transformers[10], which are in general well-suited for identifying complex
patterns. There are also papers using reinforcement learning methods such as
parameter-SARSA[11] and parameter-TD[12] in parameter learning problems.
However, we found that these methods might not be as effective in the context
of SDEs. Consider the results presented in $\bf{Figure}$ 1, for example, where
we simulated sequential data following the Heston Model $\bf{Equation}$ 33
developed in [4], and then attempted to fit it with neural networks:

Figure 1: Neural Network to Approach a Heston Model

In the upper figure we used MLP. The model is a sequential neural network with
six layers. It begins with a dense layer of 128 neurons (ReLU activation,
input size 1), followed by a 20% dropout. The next layer repeats this setup
with another 128-neuron dense layer (ReLU) and a 20% dropout. It then narrows
down to a 64-neuron dense layer (ReLU), a 32-neuron dense layer (ReLU), and
concludes with a single-neuron dense output layer. This structure is designed
for complex pattern recognition and regression tasks. The dropout layers were
included to reduce overfitting.

In the lower figure, we used another neural network, also designed for
sequential data, but starting with an input layer shaped for time series. It
features an LSTM layer with 128 units for capturing temporal dependencies,
enhanced by an attention mechanism focusing on crucial sequence parts. The
output from attention mechanism is flattened, leading to two dense layers (64
and 32 neurons, both with ReLU activation), separated by a 20% dropout layer
to reduce overfitting. The network culminates in a single-neuron output layer,
making it suitable for tasks requiring attention to specific time-dependent
patterns.

We trained the models for approximately 50 epochs, and found that the two
neural networks cannot fit the data well and are computationally very
expensive. Therefore, we sought for more efficient and effective mathematical
approaches to the problem of learning information within the realm of
stochastic differential equations. In this paper, we will introduce Kalman
filtering, a signal processing method constructed through state transition
equations, and its application to parameter estimation in various SDEs. Kalman
Filter will provide closed-form solutions for the estimation of states.

##  3 Maximum Likelihood Estimation (MLE) Framework

Parameter estimation is framed as an inverse problem, involving the deduction
of system-characterizing parameters from observational data, which is presumed
to follow the system’s rules. In scenarios where we encounter noisy data, it
is assumed that this data aligns with an equation of a predetermined form, but
the coefficients remain unknown.

###  3.1 Likelihood Function

Maximum Likelihood Estimation (MLE) is a technique for parameter estimation
when the target to estimate involves a probability distribution. In general,
let $D$ be a probability distribution with parameter $\theta$, to be
estimated. Given observations $x$ sampled from $D$, one can compute the
probability of seeing $x$ as a function of $\theta$, $P(x|\theta)$. Maximizing
this probability yields an optimal choice of $\theta$, called maximum
likelihood estimate. One could denote the function
$\mathcal{P}(\theta)=P(x|\theta)$ as the "likelihood" function of $\theta$,
although in continuous settings one need not maximize this function explicitly
in this form.

In the context of MLE applied to stochastic differential equations (SDEs), the
data given is a discrete time-series observation:
$X:\\{0,1,...,N\\}\to\mathbb{R}$ or $X=(X_{0},X_{1},...,X_{N})$. Suppose the
SDE contains an unknown parameter $\theta$. Each step in the observation
depends on the preceding step. Hence, the probability of seeing an observation
is:

$P(X|\theta)=P(X_{0}|\theta)\cdot P(X_{1}|X_{0};\theta)\cdot
P(X_{2}|X_{1};\theta)...P(X_{N}|X_{N-1};\theta)$ (1)

If $X_{t}$ for each $t$ follows a continuous conditional distribution, then
$P(X_{t}|X_{t-1};\theta)=f(X_{t}|X_{t-1};\theta)dx$ where $f$ is the
probability density function. Notice then that it suffices to drop the $dx$
and compare the product of probability density function values, instead of the
probabilities themselves. This is in turn our "likelihood function", defined
as:

$L(\theta)=f(X_{0})f(X_{1}|X_{0};\theta)\cdot
f(X_{2}|X_{1};\theta)...f(X_{N}|X_{N-1};\theta)=f(X_{0})\prod_{t=1}^{N}f(X_{t}|X_{t-1};\theta)$
(2)

The goal in an MLE, then, is to compute and maximize $L(\theta)$. In practice,
the precise formula could overflow, so instead we maximize the log-likelihood
function. Simply take the logarithm of $\bf{Equation}$ 2 to define:

$\ell(\theta)=\log(f(X_{0}))+\sum_{t=1}^{N}\log(f(X_{t}|X_{t-1};\theta))$ (3)

###  3.2 Conditional Density for Itô Diffusion

Suppose observation ${X}=(X_{0},X_{1},...,X_{N})$ follows an Itô diffusion of
form and

$dX_{t}=\mu(X_{t},t,\theta)dt+\sigma(X_{t},t,\theta)dW_{t}$ (4)

Discretizing $\bf{Equation}$ 4 gives rise to the following analogue:

$X_{t}-X_{t-1}=\mu(X_{t-1},t-1,\theta)d_{t-1}+\sigma(X_{t-1},t-1,\theta)(W_{t}-W_{t-1})$
(5)

where $dt$ is the time step. We study the conditional density
$f(X_{t}|X_{t-1};\theta)$. By definition of Brownian motion, given the
observation $X_{t-1}$, $X_{t}$ follows a normal distribution centered at
$X_{t-1}+\mu(X_{t-1},t-1,\theta)d_{t}$, with variance
$(\sigma(X_{t-1},t-1,\theta))^{2}\cdot dt$. Hence we can write down the
formula:

$\displaystyle f(X_{t}|X_{t-1};\theta)=\frac{1}{\sqrt{2\pi
dt}\sigma(X_{t-1},t-1,\theta)}\exp\left\\{-\frac{1}{2}\frac{(X_{t}-X_{t-1}-\mu(X_{t-1},t-1,\theta)dt)^{2}}{(\sigma(X_{t-1},t-1,\theta))^{2}\cdot
d_{t}}\right\\}$ (6)

To maximize $\bf{Equation}$ 2 with the conditional density given by
$\bf{Equation}$ 6 or $\bf{Equation}$ 3.

###  3.3 Marginal likelihood

The Kalman filter, akin to the Hidden Markov Model, differs in its use of
Gaussian distributed continuous variables instead of discrete states and
observations [13]. It’s important to assess the probability of the filter
producing a specific observed signal based on its parameters, such as prior
distribution, transformation, observation model, and control inputs. This
probability, known as the marginal likelihood, integrates over hidden state
variables, allowing computation with just the observed signal. It’s valuable
for evaluating parameter setups and comparing the Kalman filter with other
models using Bayesian methods.

The Kalman filter encapsulates a Markov process, wherein all pertinent
information from prior observations is encapsulated within the current state
estimate. The Marginal likelihood can be taken as the product of the
probabilities of each observation given the previous observations according to
the chain rule. If we have random observation state
$(z_{1},z_{2},...,z_{t},...)$, following the regular Kalman filter System
($H_{t}$ is the observation model, $R_{t}$ is the covariance of the
observation noise, and $P_{t}$ is the a posteriori estimate covariance
matrix), then:

$f(z)=\prod_{t=0}^{N}\int f\left(z_{t}\mid X_{t}\right)p\left(X_{t}\mid
z_{t-1},\ldots,z_{0}\right)dX_{t}=\prod_{t=0}^{N}\mathcal{N}\left(z_{t};H_{t}\hat{X}_{t\mid
t-1},R_{t}+H_{t}P_{t\mid t-1}H_{t}^{\top}\right)$ (7)

Consider $t=0$, we have a series of Gaussian densities, each linked to an
observation $z_{t}$ under the filtering distribution $H_{t}\hat{x}_{t\mid
t-1},R_{t}+H_{t}P_{t\mid t-1}H_{t}^{\top}$. A recursive update simplifies this
process. To avoid numerical underflow, the log marginal likelihood, $\ell=\log
p(z)$, with $\ell^{(-1)}=0$ , is computed instead. This is achieved through a
recursive formula, considering $d_{y}$ represents the dimension of the
measurement vector [14]:

$\ell^{(t)}=\ell^{(t-1)}-\frac{1}{2}\left(\tilde{y}_{t}^{\top}(R_{t}+H_{t}P_{t\mid
t-1}H_{t}^{\top})^{-1}\tilde{y}_{t}+\log\left|R_{t}+H_{t}P_{t\mid
t-1}H_{t}^{\top}\right|+d_{y}\log 2\pi\right)$ (8)

##  4 Kalman Filter Framework

One topic in control theory is to accurately estimate the state of a dynamical
system. Often, the state produces an observable measurement. The simplest
setup of such a system is characterized in discrete updating schemes:

$\displaystyle x_{t}$ $\displaystyle=A_{t-1}x_{t-1}+G_{t-1}W_{t-1}$ (9)
$\displaystyle y_{t}$ $\displaystyle=H_{t}x_{t}+\epsilon_{t}$

The first equation is the updating scheme for the state variable $x_{i}$, and
the second equation represents how the state variable $x_{i}$ produces a
measurement $y_{i}$. $W_{t}$ and $\epsilon_{t}$ are random variables,
representing noises involved in the state update and measurement processes. In
our studies, the two noises at the same time step could be correlated, but any
noise variables at different time steps are uncorrelated. Both $x$ and $y$ are
random variables. One usually refers to their expectations when talking about
"estimates".

Kalman filtering is a technique which uses real-time measurements to correct
state estimations. Let $Q_{t}$ be the covariance of process noise $W_{t}$ and
$R_{t}$ be the covariance of measurement noise $\epsilon_{t}$. Denote $P_{t}$
as the a posteriori state estimation error covariance, and $P_{t}^{-}$ as the
a priori state estimation error covariance. [15] and [16] describe the
algorithm to achieve optimal state estimation. We reproduce it below, $K_{t}$
is called the Kalman gain, which theoretically minimizes the a posteriori
estimate.

Algorithm 1 Kalman Filter

1:Input: The system matrices $A_{t}$, $H_{t}$, $Q_{t}$, $R_{t}$, and the
initial state $\hat{x}_{0}$ and covariance $P_{0}$.

2:Output: Updated state estimate $\hat{\mathbf{x}}_{t}$ and error covariance
$P_{t}$.

3:for each time step $t=1,2,\ldots$ do

4: Propagate the a priori state estimate:
$\hat{x}_{t}^{-}=\mathbb{E}(x_{t|t-1})=A_{t-1}\hat{x}_{t-1}$

5: Produce the estimated measurement:
$\hat{y}_{t}=\mathbb{E}(y_{t|t-1})=H_{t}\hat{x}_{t}^{-}$

6: Compute the Kalman gain matrix:
$K_{t}=P_{t}^{-}H_{t}^{T}[H_{t}P_{t}^{-}H_{t}^{T}+R_{t}]^{-1}$

7: Correct the a priori estimate with the actual measurement:
$\hat{x}_{t}=\hat{x}_{t}^{-}+K_{t}(y_{t}-\hat{y}_{t})$

8: Update error covariance matrix: $P_{t}=(I-K_{t}H_{t})P_{t}^{-}$, and
$P_{t+1}^{-}=A_{t}P_{t}A_{t}^{T}+G_{t}Q_{t}G_{t}^{T}$

9:end for

Delving into the realm of control systems, we can unveil the one-dimensional
formulation of the Kalman gain:

$\displaystyle\hat{x}_{t}$
$\displaystyle=\hat{x}_{t}^{-}+K_{t}\left(z_{t}-H_{t}\hat{x}_{t}^{-}\right)=\hat{x}_{t}^{-}+K_{t}\left(z_{t}-H_{t}x_{t}+H_{t}e_{t}^{-}\right)$
(10) $\displaystyle e_{t}$
$\displaystyle=x_{t}-\hat{x}_{t}=e_{t}^{-}-K_{t}\left(\epsilon_{t}+H_{t}e_{t}^{-}\right)$

Then, we update the error covariance matrix:

$P_{t}=E\left[e_{t}^{2}\right]=P_{t}^{-}+K_{t}^{2}\left(R_{t}+H_{t}^{2}P_{t}^{-}+0\right)-2K_{t}H_{t}P_{t}^{-}$
(11)

By differentiating with respect to $K_{t}$ and ingeniously setting this
derivative equal to zero, we arrive at the solution:

$K_{t}=\frac{P_{t}^{-}H_{t}}{H_{t}^{2}P_{t}^{-}+R_{t}}$ (12)

###  4.1 Linear Kalman Regression

In machine learning, linear regression is supervised learning, which assumes
that each measurement has the same weight. Kalman filtering incorporates the
consideration of internal changes in the system. The process equation is used
to predict the system’s state at the next moment based on its current state.
When we need to fit a known linear system, we can use linear Kalman regression
[17]. By changing the measured noise, system noise will affect the
effectiveness of Kalman filtering without affecting the linear regression
filtering, while changing the order of least squares will affect it.

Consider a linear regression defined in the following formula:

$Y=\alpha+\beta X+\epsilon$ (13)

where $Y$ and $X$ are response and predictor variables, respectively.
$\epsilon$ is the noise in system, for which $\text{Std}(\epsilon)=\sigma$ and
$\text{Cov}(\epsilon,X)=0$. In the Kalman Filter setting, the state at time
$t$ is: $x_{t}=\left(\begin{array}[]{c}x_{\alpha_{t}}\\\
x_{\beta_{t}}\end{array}\right)$. View the measurement at the next time step
as response variable $Y$, and then the system satisfies the following process
and measurement equations:

$\displaystyle\left(\begin{array}[]{c}x_{\alpha_{t}}\\\
x_{\beta_{t}}\end{array}\right)$
$\displaystyle=A_{t-1}\left(\begin{array}[]{c}x_{\alpha_{t-1}}\\\
x_{\beta_{t-1}}\end{array}\right)+G_{t-1}W_{t-1}\text{ with
}G_{t-1}W_{t-1}\sim\mathcal{N}\left(0,Q_{t-1}\right)\text{ and
}Q_{t-1}=\left(\begin{array}[]{cc}\sigma_{\alpha}^{2}&0\\\
0&\sigma_{\beta}^{2}\end{array}\right)$ (14) $\displaystyle y_{t}$
$\displaystyle=H_{t}\left(\begin{array}[]{c}x_{\alpha_{t}}\\\
x_{\beta_{t}}\end{array}\right)+\epsilon_{t}\text{ with
}H_{t}=\left(\begin{array}[]{ll}1&x_{t}\end{array}\right)\text{ and
}\epsilon_{t}\sim\mathcal{N}\left(0,\sigma_{\epsilon}^{2}\right)$

Normally, for a simple linear system, we can regard the transit matrix $A_{t}$
as an identity matrix. And the noise term
$G_{t}W_{t}=\left(\begin{array}[]{c}\nu^{\alpha}_{t}\\\
\nu^{\beta}_{t}\end{array}\right)$, where $\nu^{\alpha}_{t}$ and
$\nu^{\beta}_{t}$ are independent. Then the rest of the calculation is
supposed as the same as the regular Kalman Filter.

In the above $\bf{Equation}$ 14, the noise term $G_{t}W_{t}$ and
$\epsilon_{t}$ contain parameters $Q_{t}$ and $\sigma_{\epsilon}$. The idea is
to calibrate with the real-world data by Maximum Likelihood Estimation (MLE),
which will be introduced in Chapter 4. By following the marginal likelihood
function $\bf{Equation}$ 8:

$\log L\left(\alpha,\beta,Q_{k},\sigma_{\epsilon}\mid
y_{t},y_{t-1},\ldots,y_{1},y_{0}\right)=-\frac{1}{2}\sum_{i=1}^{N}\left(\frac{r_{t}^{2}}{S_{t}}+\log
S_{t}+\log 2\pi\right)$ (15)

where $r_{t}$ is the error between the true and estimated state. $S_{t}$ is
the variance of the system $\bf{Equation}$ 13, the system contains mean:
$\text{E}[y_{t}|\alpha_{t-1,t-1}\beta_{t-1,t-1}]=\alpha_{t,t-1}+\beta_{t,t-1}X_{t-1}+\epsilon_{t-1}$
and variance:
$\text{Var}[y_{t}|\alpha_{t-1,t-1}\beta_{t-1,t-1}]=x_{t}^{2}P_{t,t-1}+\sigma_{\epsilon}^{2}$.

###  4.2 Extended Kalman Filter (EKF)

In the above discussion, the state and measurement equations are assumed to be
linear. However, in reality, these equations are often nonlinear. The Extended
Kalman Filter (EKF) is a technique that generalizes Kalman filtering to such
settings.

Suppose the state and measurement equations for a dynamical system are:

$\displaystyle x_{t}$ $\displaystyle=f(x_{t-1},w_{t-1})\quad w_{t-1}\sim
N\left(0,Q_{t-1}\right)\hskip 10.00002pt\text{state}$ (16) $\displaystyle
y_{t}$ $\displaystyle=h(y_{t-1},\epsilon_{t-1})\quad\epsilon_{t-1}\sim
N\left(0,R_{t-1}\right)\hskip 13.00005pt\text{measurement}$

Where $f,h$ are assumed to be smooth functions. To linearize this system, we
formulate a benchmark state trajectory and study the difference between the
actual state and the benchmark state. $Q_{t}$ and $R_{t}$ is the covariance
matrix of two noise terms. In addition, $w_{t-1}$ is uncorrelated with $x_{t}$
and $\epsilon_{t-1}$ is uncorrelated with $y_{t}$.

For a function $f(x)$, we know that its first-order Taylor expansion form is:

$f(x)=f\left(x_{k}\right)+\left(x-x_{k}\right)f^{\prime}\left(x_{k}\right)$
(17)

As a state function, x estimates its linear prior process is
$\hat{x}_{t}^{-}=E\left[x_{t}\right]$, and
$\hat{x}_{t}=E\left[x_{t}|y_{t}\right]$. Assuming that the new input $u_{t}$
arrives and the true value of the state is $x_{t-1}$, a first-order Taylor
expansion is used near the estimated value $\hat{x}_{t-1}$:

$\displaystyle f\left(x_{t-1},u_{t}\right)$ $\displaystyle\approx
f\left(\hat{x}_{t-1},u_{t}\right)+f^{\prime}\left(\hat{x}_{t-1},u_{t}\right)\left(x_{t-1}-\hat{x}_{t-1}\right)$
(18) $\displaystyle f^{\prime}\left(\hat{x}_{t-1},u_{t}\right)$
$\displaystyle=\partial f(x_{t-1},u_{t})/\partial x_{t-1}$

where $f^{\prime}\left(\hat{x}_{t-1},u_{t}\right)$ is a Jacobian matrix, which
is the same as matrix $A$ in Kalman filtering we introduced at the beginning
of the section. According to the same principle, based on our initial
consideration of linear Kalman filtering in $\bf{Equation}$ 9, the
corresponding Jacobian matrix can be represented as:

$\displaystyle(A_{t-1})_{ij}$ $\displaystyle=\partial f_{i}/\partial
x_{j}\left(\hat{x}_{t-1},u_{t}\right)$ $\displaystyle(W_{t-1})_{ij}$
$\displaystyle=\partial f_{i}/\partial w_{j}\left(\hat{x}_{t-1},u_{t}\right)$
(19) $\displaystyle(H_{t})_{ij}$ $\displaystyle=\partial h_{i}/\partial
x_{j}\left(\hat{x}_{k}^{-},u_{t}\right)$ $\displaystyle(\epsilon_{t})_{ij}$
$\displaystyle=\partial
h_{i}/\partial\epsilon_{j}\left(\hat{x}_{k}^{-},0\right)$

Then, under the same foundational assumptions guiding the Kalman Filter, we
meticulously repeat the process:

Algorithm 2 Extended Kalman Filter

1:Input: Initial state estimate $\hat{\mathbf{x}}_{0}$, error covariance
$P_{0}$, system matrices $A_{t}$, $H_{t}$, $Q_{t}$, $R_{t}$, and measurement
$y_{t}$.

2:Output: Updated state estimate $\hat{\mathbf{x}}_{t}$ and error covariance
$P_{t}$.

3:for each time step $t=1,2,\ldots$ do

4: Prediction step: Propagate the a priori state estimate:
$\hat{\mathbf{x}}_{t}^{-}=\mathbf{f}(\hat{\mathbf{x}}_{t-1},0)$

5: Measurement prediction step: Produce the estimated measurement:
$\hat{y}_{t}=h(\hat{\mathbf{x}}_{t}^{-},0)$

6: Kalman gain computation: Compute the Kalman gain matrix:
$K_{t}=P_{t}^{-}H_{t}^{T}(H_{t}P_{t}^{-}H_{t}^{T}+\epsilon_{t}R_{t}\epsilon_{t}^{T})^{-1}$

7: Update step: Correct the a priori estimate using the actual measurement:
$\hat{\mathbf{x}}_{t}=\hat{\mathbf{x}}_{t}^{-}+K_{t}(y_{t}-\hat{y}_{t})$

8: Error covariance update: $P_{t}=(I-K_{t}H_{t})P_{t}^{-}$ and
$P_{t}^{-}=A_{t}P_{t-1}A_{t}^{T}+W_{t}Q_{t-1}W_{t}^{T}$

9:end for

Like what we introduced in Linear Kalman Regression. If we want to find the
parameter, we have to construct the likelihood function with marginal
likelihood function in $\bf{Equation}$ 8. By adjusting the equation with EKF,
we have:

$L_{1:T}=\sum_{t=1}^{T}\left[\ln\left(P_{t}\right)+\frac{{(y_{t}-\hat{x}_{t}^{-}})^{2}}{P_{t}}\right]$
(20)

###  4.3 Particle Extended Kalman Filter (Particle-EKF)

So far, we introduced Kalman filters and EKF based on the assumption of
Gaussian distribution. However, Gaussian distribution does not remain
unchanged in nonlinearity[18]. To address this issue, EKF focuses on mean and
covariance. However, this method cannot consistently produce precise results,
highlighting the need for more detailed solutions in such complex situations.

One method to solve this is by using a lot of sample points to show what the
distribution looks like, which is known as particle filtering[18]. The
particle filter is based on Monte Carlo simulation. The use of particle sets
to represent probabilities can be applied to any form of state space model.
The core idea is to represent the distribution of random state particles
extracted from a posterior probability, which is a sequential importance
sampling method. Details about how to accomplish the particle EKF is presented
in Algorithm 3.

Particle filters do not necessarily assume Gaussian noise. Therefore, we do
not need to maximize a Gaussian-assumed likelihood function for parameter
estimation as we would with Linear Kalman Regression or the EKF. The text[19]
explains that the total likelihood of a given state, considering all past
observations, is the product of various individual likelihoods. The
logarithmic likelihood equation simplifies the process by converting
multiplication into addition, thereby facilitating easier computation[20]:

$\ln\left(L_{1:T}\right)=\sum_{t=1}^{T}\ln\left(l_{t}\right)=\sum_{t=1}^{T}\ln\left(\int
p\left(y_{t}\mid x_{t}\right)\frac{p\left(x_{t}\mid
z_{1:t-1}\right)}{q\left(x_{t}\mid x_{t-1},y_{1:t}\right)}q\left(x_{t}\mid
x_{t-1},y_{1:t}\right)dx_{t}\right)$ (21)

The likelihood $l_{t}$ for step $t$ is approximated by a sum over $N$
simulations. Each term in the sum is a ratio of the probability of observation
$y_{t}$ given the state $x_{t}^{(i)}$, and a proposal distribution $q$. This
ratio is weighted by the probability of the state $x_{t}^{(i)}$ given the
previous state $x_{t-1}^{(i)}$. This approximation is used for sequential
importance sampling weight update in particle filtering algorithms.

$l_{t}\approx\sum_{i=1}^{N}\frac{p\left(y_{t}\mid
x_{t}^{(i)}\right)p\left(x_{t}^{(i)}\mid
x_{t-1}^{(i)}\right)}{q\left(x_{t}^{(i)}\mid x_{t-1}^{(i)},y_{1:t}\right)}$
(22)

Algorithm 3 Particle Extended Kalman Filter

1:Initialize:

2:Set $t=0$, simulation particles $N$, and choose initial state $x_{0}$ and
initial covariance $P_{0}>0$

3:for $i=1$ to $N$ do

4: Generate $x_{0}^{(i)}=x_{0}+\sqrt{P_{0}}\cdot Z^{(i)}$ $\triangleright$
$Z^{(i)}$ is a standard Gaussian random number

5: Set $P_{0}^{(i)}$ to $P_{0}$ and initial weight $w_{0}^{(i)}=1/N$

6:end for

7:

8:Algorithm:

9:while $t<N$ do

10: for each simulation index $i$ do

11: Apply Kalman Filter or EKF to get $\hat{x}_{t}^{(i)}$ from
$\hat{x}_{t-1}^{(i)}$

12: $P_{t}^{(i)}$ is the associated a posteriori error covariance matrix.

13: Generate $\tilde{x}_{t}^{(i)}=\hat{x}_{t}^{(i)}+\sqrt{P_{t}^{(i)}}\cdot
Y^{(i)}$ $\triangleright$ $Y^{(i)}$ is a standard Gaussian random number

14: end for

15: for each $i$ between 1 and $N$ do

16: Calculate the weights:

17: $w_{t}^{(i)}=w_{t-1}^{(i)}\cdot\frac{p(y_{t}|\tilde{x}_{t}^{(i)})\cdot
p(\tilde{x}_{t}^{(i)}|\hat{x}_{t-1}^{(i)})}{q(\tilde{x}_{t}^{(i)}|\hat{x}_{t-1}^{(i)},y_{t})}$
$\triangleright$ $q()$ is the normal density with mean $\hat{x}_{t}^{(i)}$ and
variance $P_{t}^{(i)}$

18: end for

19: Normalize the weights:
$\hat{w}_{t}^{(i)}=w_{t}^{(i)}/\sum_{j=1}^{N}w_{t}^{(j)}$

20: Resample the points $\tilde{x}_{t}^{(i)}$ and get $x_{t}^{(i)}$

21: Reset weights $w_{t}^{(i)}=\hat{w}_{t}^{(i)}=1/N$

22: t += 1

23:end while

##  5 SDE Simulation Results

###  5.1 Ornstein–Uhlenbeck (OU) Process

The Ornstein-Unlenbeck Process is an SDE modeling interest rates. It is
defined in the following formula:

$dX_{t}=\theta(\mu-X_{t})dt+\sigma dW_{t}$ (23)

where $W_{t}$ is a Brownian motion. Parameters $\theta,\mu,\sigma$ are
positive constants and to be estimated. Note that the $\theta$ here is not to
be confused with $\theta$ in previous sections, which represents a general
vector of parameters. The SDE could be solved analytically as:
$X_{t}=X_{0}e^{-\theta t}+\mu(1-e^{-\theta t})+\displaystyle\int_{0}^{t}\sigma
e^{\theta(s-t)}dW_{s}$

####  5.1.1 MLE of OU Process

Discretizing the SDE following $\bf{Equation}$ 5 gives:

$X_{t}=X_{t-1}+\theta(\mu-X_{t-1})dt+\sigma Z_{t-1}$

where $dt$ is the time step, and $Z_{i}$ follows independently a normal
distribution $N(0,dt)$

One computes the conditional density $f(X_{t}|X_{t-1});\theta)$ using
$\bf{Equation}$ 6 and arrives at:

$f(X_{t}|X_{t-1})=\frac{1}{\sqrt{2\pi
dt}\sigma}\exp\left\\{-\frac{1}{2}\frac{[(X_{t}-X_{t-1}-\theta(\mu-
X_{t-1})]^{2}}{\sigma^{2}\cdot dt}\right\\}$ (24)

In numerical studies, we first fix a starting point $X_{0}$ of the SDE. We
manually input parameters $\theta_{0},\mu_{0},\sigma_{0}$ and generate
simulated data. Our time step is 0.499, and we propagated the process for 1000
iterations. Then, based on the simulated data, we optimize the log-likelihood
function $\bf{Equation}$ 3 over choices of $\theta,\mu,\sigma$ using the above
calculations with scientific computing software. We aim at accurately
reconstructing the original parameters. Afterwards, we use the estimated
parameters to reconstructed another set of data. Details are available in our
code on GitHub.

Our original parameters are: $\theta=1,\mu=2,\sigma=3$. The reconstructed
parameters are: $\theta=1.006,\mu=2.065,\sigma=3.187$ with running time:
4.418s. Plotting the original and reconstructed processes yields:

Figure 2: Simulation of Ornstein–Uhlenbeck Process with original and MLE

In the figure, the blue line is the original simulated data, and the red line
is the reconstructed data. The horizontal axis is time, the vertical axis is
the value of $X_{t}$.

We also applied MLE to estimate the Black–Karasinski (B-K) model, which is
mainly used for the pricing of exotic interest rate derivatives such as
American and Bermudan bond options. The SDE of the B-K model is given by:

$d\ln(r)=(\theta-\alpha\ln(r))dt+\sigma dW_{t}$ (25)

Similarly, using $\bf{Equation}$ 6, the likelihood function of B-K model can
be written as:

$f(\ln(r_{t})|\ln(r_{t-1}))=\frac{1}{\sqrt{2\pi
dt}\sigma}\exp\left\\{-\frac{1}{2}\frac{[(\ln(r_{t})-\ln(r_{t-1})-(\theta-\alpha\ln(r_{t}))dt]^{2}}{\sigma^{2}\cdot
dt}\right\\}$ (26)

Our original parameters are set to $\theta=1.0$, $\alpha=0.8$, $\sigma=0.6$,
$dt=1$. MLE estimates the reconstructed parameters as: $\theta=0.992$,
$\alpha=0.801$, $\sigma=0.587$. Plotting the original and reconstructed
processes yields:

Figure 3: Simulation of Black–Karasinski Model with original and MLE

####  5.1.2 Kalman Filter of OU Process

Following linear system we introduced in Linear Kalman Filter, $\bf{Equation}$
13, the State and Measurement equations are:

$\displaystyle x_{t}$ $\displaystyle=\alpha+\beta
x_{t-1}+G_{t}W_{t}\quad\text{ with }\quad
G_{t}W_{t}\sim\mathcal{N}\left(0,\sigma_{GW}^{2}\right)$ (27) $\displaystyle
y_{t}$ $\displaystyle=x_{t}+\epsilon_{t}\quad\text{ with
}\quad\epsilon_{t}\sim\mathcal{N}\left(0,\sigma_{\epsilon}^{2}\right)$

By the defined formula of OU process shows above, setting: $\alpha=\theta\mu
dt,\text{ }\beta=1-\theta dt,\text{ }\text{var}[GW]=\sigma^{2}dt,\text{ and
}\text{var}[\epsilon]=1\times 10^{-4}$. And following $\bf{Equation}$ 14, the
new State and Measurement equations are:

$\displaystyle\left(\begin{array}[]{c}1\\\ x_{t}\end{array}\right)$
$\displaystyle=A_{t}\left(\begin{array}[]{c}1\\\
x_{t-1}\end{array}\right)+\left(\begin{array}[]{c}0\\\
G_{t}W_{t}\end{array}\right)\text{ with }A_{t}=\left(\begin{array}[]{cc}1&0\\\
\alpha&\beta\end{array}\right),\text{ and
}G_{t}W_{t}\sim\mathcal{N}\left(0,\sigma_{GW}^{2}\right)$ (28) $\displaystyle
y_{t}$ $\displaystyle=H_{t}\left(\begin{array}[]{c}1\\\
x_{t}\end{array}\right)+\epsilon_{t}\quad\text{ with }\quad H_{t}=1\text{ and
}\epsilon_{t}\sim\mathcal{N}\left(0,\sigma_{\epsilon}^{2}\right)$

Predict steps:

$\hat{x}_{t\mid t-1}=\alpha+\beta\hat{x}_{t-1\mid t-1}\quad\text{ and }\quad
P_{t\mid t-1}=\beta^{2}P_{t-1\mid t-1}+\sigma_{GW}^{2}\text{. }$ (29)

Supporting Variables:

$\displaystyle r_{k}=y_{t}-\hat{x}_{t\mid t-1},\quad S_{t}=P_{t\mid
t-1}+\sigma_{\epsilon}^{2},\quad K_{t}=\frac{P_{t\mid t-1}}{S_{t}}$ (30)

Update step:

$\hat{x}_{t\mid t}=\hat{x}_{t\mid t-1}+K_{t}r_{t}\quad\text{ and }\quad
P_{t\mid t}=P_{t\mid t-1}\left(1-K_{t}\right)$ (31)

In each iteration, the log likelihood function is calculated in
$\bf{Equation}$ 15.

Figure 4: OU Process Tracking with Kalman Filter

$\bf{Figure}$ 4 shows the tracking normal OU process with Kalman Filter with
same guesses parameter as MLE. The root mean square error (RMSE) between the
two states is 4.560e-05, which performs a good fitting.

By minimizing the log-likelihood function with Limited-memory BFGS, the
reconstructed parameters are: $\theta=0.9267$, $\mu=1.996$, $\sigma=2.902$
with the original parameters are as same as the previous OU Process part (
$\theta=1,\mu=2,\sigma=3$). In terms of estimating the parameters, there is no
significant difference between the two results, but the parameter estimation
method using Kalman filtering only runs for 0.30717s.

###  5.2 OU Process with Jump

For application purposes, it is often realistic to add a jump process in
addition to Brownian noise. As an example, the discretized OU process can be
adapted to

$X_{t}=X_{t-1}+\theta(\mu-X_{t-1})dt+\sigma Z_{t-1}+J_{t-1}$ (32)

where each $J_{i}$ takes parameters $\lambda_{j},\mu_{j},\sigma_{j}$ and
behaves identically and independently in the following way:

  1. 1.

sample a random number $\gamma\sim N(0,1)$

  2. 2.

if $\gamma<\lambda_{j}$, $J_{i}$ is "active", i.e., sample $J_{i}\sim
N(\mu_{j},\sigma_{j}^{2})$

  3. 3.

if $\gamma\geq\lambda_{j}$, $J_{i}=0$

We remark that generally jump processes could be much more complicated.

####  5.2.1 MLE of OU Process with Jump

With the above example, one can in fact still analytically compute the log-
likelihood function. This results from the design of $J_{i}$ by our
assumption.

Following our previously introduced framework, we focus on computing the
conditional density $f(X_{t}|X_{t-1})$. This is done by considering the two
cases described in bullet points 2 and 3 above. If $J_{t-1}$ is active (case
2), $X_{t}$ follows a normal distribution centered at
$\left(X_{t-1}+\theta(\mu-X_{t-1})dt+\mu_{j}\right)$, with variance
$\sigma^{2}dt+\sigma_{j}^{2}$. If $J_{t-1}$ is identically zero (case 3),
$X_{t}\sim N\left(X_{t-1}+\theta(\mu-X_{t-1})dt,\sigma^{2}dt\right)$ as it is
in the simple OU process.

With this established, knowing the density for normal distributions, we
compute:

$\displaystyle P(X_{t}|X_{t-1})=$
$\displaystyle=P(X_{i+1},J_{i}=0|X_{i})+P(X_{i+1},J_{i}\neq 0|X_{i})$
$\displaystyle=P(X_{i+1}|J_{i}=0,X_{i})\cdot P(J_{i}=0)+P(X_{i+1}|J_{i}\neq
0,X_{i})\cdot P(J_{i}\neq 0)$
$\displaystyle=(1-F_{N}(\lambda_{j}dt))\cdot\frac{1}{\sqrt{2\pi
dt}\sigma}\cdot\exp\left\\{-\frac{1}{2}\left(\frac{X_{i+1}-[X_{i}+\theta(\mu-
X_{i})dt]}{\sigma\sqrt{dt}}\right)^{2}\right\\}dx$
$\displaystyle+F_{N}(\lambda_{j}dt)\cdot\frac{1}{\sqrt{2\pi(\sigma^{2}dt+\sigma_{j}^{2})}}\cdot\exp\left\\{-\frac{1}{2}\left(\frac{X_{i+1}-[X_{i}+\theta(\mu-
X_{i})dt+\mu_{j}]}{\sqrt{\sigma^{2}dt+\sigma_{j}^{2}}}\right)^{2}\right\\}dx$

where we used Bayes’ Theorem in the second step. $F_{N}$ denotes the
cumulative distribution function for the standard normal distribution, and the
$dx$ is the technicality term for calculating probability. Omitting the $dx$
gives the conditional density $f(X_{t}|X_{t-1})$.

With such, we repeated the procedure in $\bf{Section}$ 5.1.1 to examine the
effectiveness of MLE in this setting. We enforced a bound of $(10^{-15},6)$
for all parameters during optimization. The results are summarized in the
following with running time: 27.83828s:

Parameter | $\theta$ | $\mu$ | $\sigma$ | $\lambda_{j}$ | $\mu_{j}$ | $\sigma_{j}$
---|---|---|---|---|---|---
Original | 1 | 2 | 4 | 0.5 | 1 | 1
Reconstructed | 0.86 | 0 | 4.005 | 0 | 1.221 | 0.798
Table 1: MLE results for OU process with jump

In the following Figure 5, the blue lines represent original simulated data,
and the red lines are reconstructed data.

Figure 5: Simulation of Ornstein–Uhlenbeck Process with jump and MLE

####  5.2.2 Kalman Filter of OU Process with Jump

In the example of the OU process with a jump, given its linear nature, the
Linear Kalman Filter can be effectively employed to describe the state and
measurement, as detailed in $\bf{Equation}$ 27. The thing we need to aware of
that adding jump term to the noise, which means:
$\text{var}[GW]=\sigma_{GW}^{2}dt+\lambda_{j}\cdot\mu_{j}^{2}dt$. The
remaining parts of the equation, including the prediction and update steps,
follow the methodology outlined in $\bf{Section}$ 5.1.2.

Figure 6:  OU Process with Jump Tracking with Kalman Filter

The illustration in $\bf{Figure}$ 6 depicts the tracking of an OU process that
includes a jump term, utilizing a Kalman Filter with parameters identical to
those used in the MLE. The RMSE between the two states is 0.0001021,
indicating a high degree of accuracy in fitting.

By minimizing the log-likelihood function using Limited-memory BFGS, the
reconstructed parameters are obtained as:

Parameter | $\theta$ | $\mu$ | $\sigma$ | $\lambda_{j}$ | $\mu_{j}$ | $\sigma_{j}$
---|---|---|---|---|---|---
Original | 1 | 2 | 4 | 0.5 | 1 | 1
Reconstructed | 0.7858 | 6.0 | 4.844 | 0.5525 | 1.496 | 1.0
Table 2: Kalman filter results for OU process with jump

We found the results of the $\mu$ term to be inaccurate. This could be
attributed to the filters’ comparative lack of sensitivity towards the
parameter set, or perhaps due to the non-uniqueness of the optimal parameter
configuration.

The method of parameter estimation by using the solution of the Kalman filter
does not show better results very intuitively compared to the MLE. However,
compared to the 27.83828s spent by MLE, the running time of Kalman filter is
0.265414s. There is a significant difference between them.

Next, we aim to evaluate a more comprehensive model within our Kalman Filter
framework to ascertain its efficacy in more complex scenarios.

###  5.3 Heston Model

In 1993, Heston enhanced the traditional Black-Scholes (BS) model by
introducing stochastic volatility[4], thereby deriving a model tailored for
European option pricing. Heston formulated the SDE under the Physical Measure
as:

$\begin{gathered}dS_{t}=\mu S_{t}dt+\sqrt{v_{t}}S_{t}dW_{t}^{1},\\\
dv_{t}=\kappa\left(\theta-v_{t}\right)dt+\xi\sqrt{v_{t}}dW_{t}^{2},\\\
dW_{1}\cdot dW_{2}=\rho dt\end{gathered}$ (33)

In $\bf{Equation}$ 33, $S_{t}$ is the price of the asset, $v_{t}$ is the
instantaneous variance from the CIR process. $W_{t}^{1}$ and $W_{t}^{2}$ are
two wiener process with instantaneous correlation $\rho$. The parameter $\mu$
is the risk-neutral rate of return, $\kappa$ is the adjustment rate, $\theta$
is the long-term run average variance, and $\xi$ is the volatility
fluctuations [21].

By expressing $dS_{t}$ as $d\ln(S_{t})$, this form simplifies the calculation
of volatility, which typically impacts the rate of return rather than the
absolute price. After transforming by Ito’s Lemma, we get:

$d\ln\left(S_{t}\right)=\left(\mu-\frac{1}{2}v_{t}\right)dt+\sqrt{v_{t}}dW_{t}^{1}$
(34)

Since volatility cannot be negative, we use the Feller condition:
$\frac{2\kappa\theta}{\xi^{2}}>1$ [22]. Also, since pricing needs to be
guaranteed to be arbitrage-free, we need to convert the SDE under Physical
Measure to Martingale Measure, using the market price of risk:

$\displaystyle
d\ln\left(S_{t}\right)=\left(\mu_{s}-\frac{1}{2}v_{t}\right)dt+\sqrt{v_{t}}\widetilde{W}_{t}^{1},\qquad\widetilde{W}_{t}^{1}$
$\displaystyle=dW_{t}^{1}+\frac{\mu-\mu_{s}}{\sqrt{v_{t}}}dt$ (35)
$\displaystyle dv_{t}=[\kappa\left(\theta-v_{t}\right)-\lambda
v(t)]dt+\xi\sqrt{v_{t}}\widetilde{W}_{t}^{2},\qquad\widetilde{W}_{t}^{2}$
$\displaystyle=dW_{t}^{2}+\frac{\left(W_{t}^{2}-\rho
W_{t}^{1}\right)}{\sqrt{1-\rho^{2}}}dt$

then:

$v_{t}=v_{t-1}+\left(\rho\xi-\theta
v_{t-1}\right)dt+\xi\sqrt{v_{t-1}}\sqrt{dt}W_{t-1}^{2}-\rho\xi\left[\ln
S_{t-1}+\left(\mu_{s}-\frac{1}{2}v_{t-1}\right)dt+\sqrt{v_{t-1}}\sqrt{dt}W_{t-1}^{1}-\ln
S_{t}\right]$ (36)

Considering the uncorrelated with $W_{t}^{1}$ and $W_{t}^{2}$, we have the
observe and measurement for Heston:

$\displaystyle x_{t}$
$\displaystyle=v_{t}=v_{t-1}+\left[\left(\kappa\left(\theta-
v_{t}\right)-\rho\xi\mu_{s}\right)-\left(\theta-\frac{\rho\xi}{2}\right)v_{t-1}\right]dt+\rho\xi\ln\left(\frac{S_{t}}{S_{t-1}}\right)+\xi\sqrt{1-\rho^{2}}\sqrt{v_{t-1}}\sqrt{dt}\tilde{W}_{t-1}^{2}$
(37) $\displaystyle y_{t}$ $\displaystyle=\ln S_{t}=\ln
S_{t-1}+\left(\mu-\frac{1}{2}v_{t-1}\right)dt+\sqrt{v_{t-1}}\sqrt{dt}W_{t-1}^{1}\
$

Overall, the Heston model is a two-factor nonlinear model. We will fit it by
using filtering. After that, based on our filtering result, we will try to
find the parameters based on minimizing the log-likelihood function.

Figure 7:  Heston with Euler–Maruyama Method

####  5.3.1 Extended Kalman Filter in Heston

In terms of the EKF, the following state equation is presented in Heston
model:

$x_{t}=f\left(x_{t-1},w_{t-1}\right)=\left(\begin{array}[]{c}\ln
S_{t-1}+\left(\mu_{S}-\frac{1}{2}v_{t-1}\right)dt+\sqrt{v_{t-1}}\sqrt{dt}W_{t-1}^{1}\\\
v_{t-1}+\left(\kappa\left(\theta-v_{t}\right)-\theta
v_{t-1}\right)dt+\xi\sqrt{v_{t-1}}\sqrt{dt}W_{t-1}^{2}\end{array}\right)$ (38)

Following the procedure outlined shown in the input part of Algorithm 2, the
pertinent matrix is displayed below:

$Q_{t}=\left(\begin{array}[]{ll}1&\rho\\\ \rho&1\end{array}\right),\quad
A_{t}=\left(\begin{array}[]{cc}1&-\frac{1}{2}dt\\\ 0&1-\theta
dt\end{array}\right),\quad
W_{t}=\left(\begin{array}[]{cc}\sqrt{v_{t-1}}\sqrt{dt}&0\\\
0&\xi\sqrt{v_{t-1}}\sqrt{dt}\end{array}\right),\quad
w_{t}=\left(\begin{array}[]{c}W_{t}^{1}\\\ W_{t}^{2}\end{array}\right)$ (39)

We follow the formulation by Alireza Javaheri[23] and write state-space into
one dimension:

$A_{t}=1-\left(\theta-\frac{1}{2}\rho\xi\right)dt,\quad
W_{t}=\xi\sqrt{1-\rho^{2}}\sqrt{v_{k-1}}\sqrt{dt},\quad
H_{t}=-\frac{1}{2}dt,\quad Q_{t}=\sqrt{v_{t}}\sqrt{dt}$ (40)

Figure 8:  Heston Tracking with EKF

The illustration in Figure 8 depicts the tracking of a Heston model using an
EKF with parameters identical to those in the likelihood function. The RMSE
between the two states is 1.268415 and running time 0.01886s, indicating a
high degree of accuracy. Regarding parameter estimation, by minimizing the
log-likelihood function as shown in Equation 21, we obtained the estimated
parameters presented in Table 3.

| $\boldsymbol{\mu}$ | $\boldsymbol{\theta}$ | $\boldsymbol{\rho}$ | $\boldsymbol{K}$ | $\boldsymbol{\xi}$
---|---|---|---|---|---
Actual | 0.05 | 1.5 | 0.04 | 0.3 | -0.6
Guess | 0.1 | 1.0 | 0 | 0.1 | 0
Estimate | 1 | 2 | 0 | 0.5 | -1
Table 3: EKF Estimate Heston’s parameters

We found the results for the $\mu$ term to be inaccurate, though other
predictions appear reasonable. Reasonable estimates of parameters and low RMSE
values suggest that the application of the EKF in the Heston model has been
successful.

####  5.3.2 Particle Extended Kalman Filter in Heston

In Section 4.3, we introduced the Particle Extended Kalman Filter. In this
section, we adapt it to a one-dimensional state space with the following
formulation:

$\displaystyle n(x,m,s)$
$\displaystyle=\frac{1}{\sqrt{2\pi}s}\exp\left(-\frac{(x-m)^{2}}{2s^{2}}\right),\
q\left(\tilde{v}t^{(i)}\mid x{t-1}^{(i)},\ln S_{0:t-1}\right)$
$\displaystyle=n\left(\tilde{x}_{t}^{(i)},m=\hat{x}_{t}^{(i)},s=\sqrt{P_{t}^{(i)}}\right),$
(41)

where the normal density is characterized by mean $m$ and standard deviation
$s$. Consequently, we derive:

$p\left(\ln S_{t-1}\mid\tilde{x}k^{(i)}\right)=n\left(\ln
S{t-1},m=\operatorname{Ln}S_{t}+\left(\mu_{S}-\frac{1}{2}\tilde{x}_{t}^{(i)}\right)dt,s=\sqrt{\tilde{x}_{t}^{(i)}}\sqrt{dt}\right),$
(42)

and:

$\displaystyle p\left(\tilde{x}t^{(i)}\mid
x{t-1}^{(i)}\right)=n\left(\tilde{x}t^{(i)},m=x{t-1}^{(i)}+\left(\kappa\left(\theta-
x_{t-1}\right)-\rho\xi\mu_{S}\right)dt-\frac{1}{2}\rho\xi
x_{t-1}^{(i)}dt+\rho\xi\left(\ln S_{t}-\ln S_{t-1}\right),s\right).$ (43)

This provides us with the necessary densities for the filter implementation:

$s=\xi\sqrt{1-\rho^{2}}\sqrt{x_{t-1}^{(i)}}\sqrt{dt}$ (44)

With 1000 particles being used, the estimation of the observed state is given
by:

$\ln S_{t-1}^{-}=\frac{1}{1000}\sum_{i=1}^{1000}\ln\hat{S}_{t-1}^{(i)}$ (45)

obtained from the EKF applied to $x_{t-1}$. See $\bf{Algorithm3}$ for more
details.

Figure 9:  Heston Tracking with Particle-EKF

Figure 9, illustrates the tracking performance of a Heston model using the
Particle-EKF. With an RMSE score of 8.9011336, running time of 130.56615s, the
results indicate a commendable level of fit.

###  5.4 Bates Model

The bates model is based on the Heston model, which extends the volatility by
incorporating a Merton log-normal jump[24]. As $\bf{Equation}$ 46 shows below,
consequently, the number of jumps occurring in the interval between $t+dt$
will be denoted as $dN_{t}$, which is a Poisson counter.

$dS_{t}=\left(\mu_{s}+\lambda
j\right)S_{t}dt+\sqrt{v_{t}}S_{t}dW_{t}-S_{t}jdN_{t}$ (46)

Here, $\lambda$ is intensity, and j is the Gaussian random variable that has a
fixed fractional jump size $[0,1)$, which has a close definition in
$\bf{Section}$ 5.2. And by applying Ito’s lemma for semi-Martingales:

$d\ln S_{t}=\left(\mu_{s}-\frac{1}{2}v_{t}+\lambda
j\right)dt+\sqrt{v_{t}}dW_{t}+\ln(1-j)dN_{t}$ (47)

Figure 10:  Bates with Euler–Maruyama Method

As same as what we have down in Heston, we will fit this model by using
filtering. Then, we have the transition matrix:

$x_{t}=\left(\begin{array}[]{c}\delta_{0}(0)e^{-\lambda
dt}+\delta_{0}(\ln(1-j))\left(1-e^{-\lambda dt}\right)\\\
v_{k-1}+\left[\left(\kappa\left(\theta-
v_{t}\right)-\rho\xi\left(\mu_{s}+\lambda
j\right)-\left(\theta-\frac{1}{2}\rho\xi\right)v_{t-1}\right]dt+\rho\xi\left[\ln\left(\frac{S_{t}}{S_{t-1}}\right)-\mu_{t-1}\right]+\xi\sqrt{1-\rho^{2}}\sqrt{v_{t-1}}\sqrt{dt}\tilde{W}_{t-1}^{2}\right.\end{array}\right)$
(48)

####  5.4.1 Extended Kalman Filter in Bates

We have imployed a similar method in $\bf{Section5.2.2}$, with adapting EKF to
processes with jumps, by adding $\lambda j$ to all $\mu_{s}$. Following the
$\bf{Algorithm2}$, we got:

Figure 11:  Bates Tracking with EKF

The illustration in Figure 8 depicts the tracking of a Heston model using an
EKF with parameters identical to those in the likelihood function. The RMSE
between the two states is 6.37269 and running time is 0.0081274, which is not
a bad solution. Regarding parameter estimation, by minimizing the log-
likelihood function as shown in Equation 20, we obtained the estimated
parameters presented in Table 4.

| $\boldsymbol{\mu}$ | $\boldsymbol{\theta}$ | $\boldsymbol{\rho}$ | $\boldsymbol{K}$ | $\boldsymbol{\xi}$ | $\lambda$ | j
---|---|---|---|---|---|---|---
Actual | 0.05 | 1.5 | 0.04 | 0.3 | -0.6 | 10 | 0.1
Guess | 0.1 | 1.0 | 0 | 0.1 | 0 | 8 | 0.1
Estimate | 1 | 2 | 0 | 0.5 | -1 | 6 | 0.1
Table 4: EKF Estimate Bates’ parameters

However, parameter estimation results yielded by EKF are not ideal. This could
be attributed to the filters’ comparative lack of sensitivity towards the
parameter set, or perhaps due to the non-uniqueness of the optimal parameter
configuration.

####  5.4.2 Particle Extended Kalman Filter in Bates

The approach of using Particle-EKF in the Bates model is similar to our
implementation in the Heston model. The adjustment required is the addition of
$\lambda j$ to $\mu_{s}$. For instance, refer to $\bf{Equation}$42 and
$\bf{Equation}$43:

$\displaystyle p\left(\ln S_{t-1}\mid\tilde{x}k^{(i)}\right)$
$\displaystyle=n\left(\ln
S{t-1},m=\operatorname{Ln}S_{t}+\left((\mu_{s}+\lambda
j)-\frac{1}{2}\tilde{x}_{t}^{(i)}\right)dt,s=\sqrt{\tilde{x}_{t}^{(i)}}\sqrt{dt}\right),$
(49) $\displaystyle p\left(\tilde{x}t^{(i)}\mid x{t-1}^{(i)}\right)$
$\displaystyle=n\left(\tilde{x}t^{(i)},m=x{t-1}^{(i)}+\left(\kappa\left(\theta-
x_{t-1}\right)-\rho\xi(\mu_{s}+\lambda j)\right)dt-\frac{1}{2}\rho\xi
x_{t-1}^{(i)}dt+\rho\xi\left(\ln S_{t}-\ln S_{t-1}\right),s\right).$

By simplifying these elements in one dimension form, we obtain:

$\displaystyle p\left(lnS_{t+1}\mid\tilde{x}_{k}^{(i)}\right)$
$\displaystyle=n\left(lnS_{t+1},lnS_{t}+\left(\mu_{s}+\lambda
j-\frac{1}{2}\tilde{x}_{t}^{(i)}\right)dt+\tilde{\mu}_{t}^{(i)},\sqrt{\tilde{x}_{t}^{(i)}dt}\right)$
(50) $\displaystyle p\left(\tilde{x}_{t}^{(i)}\mid x_{t-1}^{(i)}\right)$
$\displaystyle=n\left(\tilde{x}_{t}^{(i)},m,s=\xi\sqrt{1-\rho^{2}}\sqrt{x_{t-1}^{(i)}}\sqrt{dt}\right)p\left(\tilde{\mu}_{t}^{(i)}\mid\mu_{t-1}^{(i)}\right)$

And then, we have the mean and particles:

$\displaystyle{m}$ $\displaystyle=x_{t-1}^{(i)}+\left[\kappa\left(\theta-
x_{t}\right)-\rho\xi\left(\mu_{s}+\lambda
j\right)-\left(\theta-\frac{1}{2}\rho\xi\right)x_{t-1}^{(i)}\right]dt+\rho\xi\ln\left(\frac{S_{k}}{S_{k-1}}\right)-\rho\xi\mu_{t-1}^{(i)}$
(51) $\displaystyle q\left(\tilde{x}_{t}^{(i)}\mid
x_{t-1}^{(i)},lnS_{1:t+1}\right)$
$\displaystyle=n\left(\tilde{x}_{t}^{(i)},\hat{x}_{t}^{(i)},P_{t}^{(i)}\right)p\left(\tilde{\mu}_{t}^{(i)}\mid\mu_{t-1}^{(i)}\right)$

And following $\bf{Algorithm3}$, we got:

Figure 12:  Bates Tracking with Particle-EKF

The figure referenced above, Figure 12, illustrates the tracking performance
of a Bates model using Particle-EKF. With an RMSE score of 10.48056 and
running time 35.559821, the results indicate a commendable level of fit.

##  6 Conclusion & Future Work

In this chapter, we will conclude with the properties of Kalman filtering and
gain inspiration through extensive exploration.

Firstly, our discussion will pivot towards the concept of sequential
processing in Kalman filtering. The characteristic of this method is to adopt
an incremental data processing strategy, processing one observation at a time.
This progressive technology not only has high computational efficiency but
also greatly benefits in real-time application and processing of large
datasets. Therefore, Kalman filtering has shown commendable performance in
fitting SDE systems. Evident from the RMSE values observed in various models,
such as the OU process, the OU process with jumps, the Heston model, and the
Heston model with jumps, Kalman filtering shows impressive fitting results.
Specifically, the RMSE values are 4.560e-05 and 0.0001021 in the OU process,
and 1.268415 and 6.37269 in the Heston models, respectively. Intriguingly,
Kalman filtering maintains its robust performance even in the presence of
jumps in Lévy processes. This characteristic of Kalman filtering offers
valuable insights for future research, particularly in addressing noise issues
in Lévy spaces.

In addition, we consider closed-form estimation in linear systems. Here, the
Kalman filter shines brightly, providing a direct closed-form solution for
state estimation. This direct method greatly reduces computational complexity,
especially compared to resource-intensive iterative optimization methods. In
addition, in terms of computational load, Kalman filters operating in linear
systems typically have a lighter burden and are linearly proportional to the
amount of data. This is very evident in comparison to MLE. The results of
parameter estimation using Kalman filtering and MLE estimation are similar for
the same SDE. But the time required for the two methods is completely
different. Kalman filter takes 0.30717 seconds and 0.265414 seconds
respectively, but MLE is fertilized for 4.418 seconds and 27.83828 seconds.
This may be because MLE typically requires more computationally intensive
iterative optimization techniques, especially in complex or large-scale
models. In summary, the Kalman filtering method is faster in parameter
estimation of SDE systems, mainly due to its sequential updating mechanism and
the existence of closed-form solutions in some cases. At the same time, MLE
requires iterations and may require larger computational complexity.

After that, we discussed the inherent assumptions and constraints. The
assumptions of the standard Kalman filter are linear and Gaussian noise. To
adapt to nonlinear models, alternative solutions such as EKF or Particle-EKF
have begun to take effect. Although these variants may increase computational
complexity, they generally maintain their effectiveness and applicability.
Meanwhile, we evaluated the robustness of the Kalman filter in SDE systems,
and conducted tests on the EKF within the Heston model using various
parameters:

Parameter | $\mu$ | $\kappa$ | $\theta$ | $\sigma$ | $\rho$ | RMSE
---|---|---|---|---|---|---
task 1 | 0.04 | 0.3 | 1.5 | -0.6 | 0.04 | 1.268415
task 2 | 0.1 | 1.0 | 0.02 | 0.1 | -0.8 | 1.01548
task 3 | 0.3 | 2.0 | 0.01 | 0.6 | -0.1 | 0.77223
task 4 | 0.23 | 0.01 | 3 | 0.2 | 0 | 1.66183
Table 5: EKF In Heston with Different Parameter

Based on the above $\bf{Table}$ 5, We find that all the values of RMSE do not
vary outrageously, and almost all of them vary around 1. This suggests that
our EKF is stable in the Heston model.

In the next stage, our attention will shift towards a more complex framework,
using Chen’s model[25] as an example. This model is known for its subtle
description of interest rate dynamics, providing a three-factor model rather
than the two-factor structure of the Heston model we have explored so far.
Based on this progress, we plan to rigorously evaluate the adaptability of
Kalman filters by integrating jump components, similar to the methods used in
the Bates model. This will enable us to more accurately evaluate the
effectiveness of Kalman filters in fitting complex financial models.

After this, our research will focus more on prediction. Our main goal is to
capitalize on the potential of sliding window techniques as a tool for
prediction. We aim to develop a more nuanced understanding of predictive
dynamics. This forecasting journey will allow us to rigorously apply the
Kalman Filter to the S&P 500 in particular. Our intention is not just to apply
this technique in a theoretical vacuum, but to utilize it in the challenging
and unpredictable realm of real-world financial data. This will allow us to
carefully evaluate the robustness and predictive power of the filter in an
environment filled with complexities. As we progress, our research will extend
to evaluating the adaptability and accuracy of Kalman filters under a variety
of market conditions, including periods of high volatility and market
corrections. We will investigate how the filter responds to sudden economic
events and changes in market trends. In addition, we aim to compare the
performance of Kalman filters with other forecasting models and assess their
relative strengths and weaknesses in predicting market volatility.

## Github Link

All code pertaining to this project can be accessed at:

https://github.com/bao021600/Application-of-Kalman-Filter-in-Stochastic-
Differentional-Equations

## References

  * [1] BENI LAUTERBACH and PAUL SCHULTZ.  Pricing warrants: An empirical study of the black-scholes model and its alternatives.  The Journal of Finance, 45(4):1181–1209, September 1990.
  * [2] Richard Jordan, David Kinderlehrer, and Felix Otto.  The variational formulation of the fokker–planck equation.  SIAM Journal on Mathematical Analysis, 29(1):1–17, 1998.
  * [3] Nguyen Huu Du and Vu Hai Sam.  Dynamics of a stochastic lotka–volterra model perturbed by white noise.  Journal of Mathematical Analysis and Applications, 324(1):82–97, 2006.
  * [4] Steven L. Heston.  A closed-form solution for options with stochastic volatility with applications to bond and currency options.  The Review of Financial Studies, 6, no. 2:327–343, 1993.
  * [5] Ximei Wang, Xingkang He, Yanlong Zhao, and Zhiqiang Zuo.  Parameter estimations of heston model based on consistent extended kalman filter.  IFAC-PapersOnLine, 50(1):14100–14105, 2017.  20th IFAC World Congress.
  * [6] T. Kailath.  A view of three decades of linear filtering theory.  IEEE Transactions on Information Theory, 20(2):146–181, 1974.
  * [7] Haixia Kang, Yuping Wu, Lan Chen, and Xuelian Zhang.  Research on device modeling technique based on mlp neural network for model parameter extraction.  Applied Sciences, 12(3), 2022.
  * [8] K.E. ArunKumar, Dinesh V. Kalaga, Ch. Mohan Sai Kumar, Masahiro Kawaji, and Timothy M. Brenza.  Comparative analysis of gated recurrent units (gru), long short-term memory (lstm) cells, autoregressive integrated moving average (arima), seasonal autoregressive integrated moving average (sarima) for forecasting covid-19 trends.  Alexandria Engineering Journal, 61(10):7585–7603, 2022.
  * [9] Wei Song, Chao Gao, Yue Zhao, and Yandong Zhao.  A time series data filling method based on lstm—taking the stem moisture as an example.  Sensors, 20(18):5045, September 2020.
  * [10] Mohamed F. Kotb, Attia A. El-Fergany, and Eid A. Gouda.  Estimation of electrical transformer parameters with reference to saturation behavior using artificial hummingbird optimizer.  Scientific Reports, 12(1), November 2022.
  * [11] Lingwei Zhu, Go Takami, Mizuo Kawahara, Hiroaki Kanokogi, and Takamitsu Matsubara.  Alleviating parameter-tuning burden in reinforcement learning for large-scale process control.  Computers & Chemical Engineering, 158:107658, 2022.
  * [12] Steven J. Bradtke and Andrew G. Barto.  Machine Learning, 22(1/2/3):33–57, 1996.
  * [13] L.R. Rabiner.  A tutorial on hidden markov models and selected applications in speech recognition.  Proceedings of the IEEE, 77(2):257–286, 1989.
  * [14] Helmut Lütkepohl.  Introduction to Multiple Time Series Analysis.  Heidelberg: Springer-Verlag Berlin. p. 435, 1991.
  * [15] Harry Asada.  Lecture_5.pdf: Identification, estimation, and learning: Mechanical engineering, 2006.
  * [16] Harry Asada.  Lecture_6.pdf: Identification, estimation, and learning: Mechanical engineering, 2006.
  * [17] H.W. SORENSON.  Kalman filtering techniques.  volume 3 of Advances in Control Systems, pages 219–292. Elsevier, 1966.
  * [18] Katalin György, András Kelemen, and László Dávid.  Unscented kalman filters and particle filter methods for nonlinear state estimation.  Procedia Technology, 12:65–74, 2014.  The 7th International Conference Interdisciplinarity in Engineering, INTER-ENG 2013, 10-11 October 2013, Petru Maior University of Tirgu Mures, Romania.
  * [19] ALIREZA JAVAHERI.  Inside Volatility Arbitrage.  John Wiley & Sons,pages 99–102, 2005.
  * [20] Surya T. Tokdar and Robert E. Kass.  Importance sampling: a review.  WIREs Computational Statistics, 2(1):54–60, December 2009.
  * [21] Milan Mrázek and Jan Pospíšil.  Calibration and simulation of heston model.  Open Mathematics, 15(1):679–704, 2017.
  * [22] Bernt K. Øksendal.  Stochastic Differential Equations: An Introduction with Applications (Sixth ed.). Berlin: Springer. ISBN 3-540-04758-1., 2003.
  * [23] ALIREZA JAVAHERI.  Inside Volatility Arbitrage.  John Wiley & Sons,pages 87–89, 2005.
  * [24] David S. Bates.  Jumps and Stochastic volatility: Exchange Rate Processes Implicity in Deutsche Mark Options.  The Review of Financial Studies, volume 9, number 1,pages 69–107, 1996.
  * [25] Hsing-Hui Chu, Tai-Liang Chen, Ching-Hsue Cheng, and Chen-Chi Huang.  Fuzzy dual-factor time-series for stock index forecasting.  Expert Systems with Applications, 36(1):165–171, 2009.