# Infinite Nature:Perpetual View Generation of Natural Scenes from a Single
Image

Andrew Liu Richard Tucker††footnotemark:  Varun Jampani
Ameesh Makadia Noah Snavely Angjoo Kanazawa
Google Research Equal Contribution

###### Abstract

We introduce the problem of perpetual view generation—long-range generation of
novel views corresponding to an arbitrarily long camera trajectory given a
single image. This is a challenging problem that goes far beyond the
capabilities of current view synthesis methods, which work for a limited range
of viewpoints and quickly degenerate when presented with a large camera
motion. Methods designed for video generation also have limited ability to
produce long video sequences and are often agnostic to scene geometry. We take
a hybrid approach that integrates both geometry and image synthesis in an
iterative ‘ _render_ , _refine_ and _repeat_ ’ framework, allowing for long-
range generation that cover large distances over hundreds of frames. Our
approach can be trained from a set of monocular video sequences without any
manual annotation. We propose a dataset of aerial footage of natural coastal
scenes, and compare our method with recent view synthesis and conditional
video generation baselines, showing that it can generate plausible scenes for
much longer time horizons over camera trajectories covering a large distance
compared to existing methods. Please visit our project page at
https://infinite-nature.github.io/.

Figure 1: Perpetual View Generation. From a collection of aerial videos of
nature scenes (left), we propose a method that can perpetually generate novel
views for a camera trajectory covering a long distance from a single image
(right). Our method can successfully generate hundreds of frames of an aerial
video from a single input image (up to 500 shown here).  ††* indicates equal
contribution

##  1 Introduction

Consider the input image of a coastline in Fig. Infinite Nature:Perpetual View
Generation of Natural Scenes from a Single Image. Imagine flying through this
scene on a small airplane. Initially, we would see objects grow in our field
of view as we approach them. Beyond, we might find a wide ocean or new
islands. At the shore, we might see cliffs or beaches, while inland there
could be mountains or forests. As humans, we are good at imagining a plausible
world from a single picture, based on our own experiences. How can we build a
system that can similarly imagine the world beyond the edges of a single
image?

We introduce the problem of perpetual view generation: given a single image of
a scene, the goal is to continually generate new views of the scene
corresponding to an arbitrarily long camera trajectory, giving the effect of
flying through an imaginary world generated from the given image. Solving this
problem requires extrapolating new content for unseen regions and synthesizing
new details in existing regions as the camera approaches them. Building an
infinitely generative scene model has applications in content creation, novel
photo interactions, and methods that use learned world models like model-based
reinforcement learning.

However, generating a long video from a static image is an extremely
challenging problem. Two active areas of research, video synthesis and view
synthesis, both fail to scale to this problem for different reasons. Recent
video synthesis methods generate limited numbers of novel frames (e.g., 25
[40] or 48 frames [9]), even when trained with massive compute resources.
These methods take advantage of recent developments in image synthesis [20]
and apply them to the temporal domain or rely on recurrent models [10].
However, they often neglect an important element of the video’s structure—a
video is a function of both the underlying scene _and_ camera geometry. Proper
geometry imposes constraints on how a video can evolve and is crucial for
synthesizing moving camera sequences such as the ones we tackle.

In contrast, many view synthesis methods do take advantage of geometry to
synthesize high-quality novel views. However, these approaches can only
operate within a limited range of camera motions. As shown in Fig. 6, if the
camera strays too far, such methods fail catastrophically. To successfully
generate distant views, an algorithm will need to inpaint disoccluded regions,
extrapolate (i.e. outpaint) unseen regions beyond the boundary of the previous
frames, and add detail to (super-resolve) areas that approach the camera over
time.

We propose a hybrid framework that takes advantage of both geometry and image
synthesis techniques to address these challenges. Specifically, we use
disparity maps to encode the scene geometry, and decompose the perpetual view
generation task into the framework of _render_ -_refine_ -and-_repeat_. First,
we render the current frame from a new viewpoint, using disparity to ensure
that scene content moves in a geometrically correct manner. Then, we refine
the resulting image and geometry. This step adds detail and synthesizes new
content in areas that require inpainting, outpainting, and super-resolution.
Because we refine both the image and disparity, the whole process can be
repeated in an auto-regressive manner, allowing for perpetual generation of
novel views.

To train our system, we curated a large dataset of drone footage of nature and
coastal scenes from over 700 videos, spanning 2 million frames. We run a
structure from motion pipeline to recover 3D camera trajectories, and refer to
this as the Aerial Coastline Imagery Dataset (ACID), which we have released
publicly. Our trained model can generate sequences of hundreds of frames while
maintaining the aesthetic feel of an aerial coastal video, even though after
just a few frames, the camera has moved beyond the edges of the original
image.

Our experiments show that the render-refine-repeat structure of our framework
is key to tackling this problem. Compared to recent view synthesis and video
generation baselines, we show that our approach can produce plausible frames
for much longer time horizons than prior methods. While there is more to be
done on this immensely challenging task, our work shows the potential of
merging geometry and generative models, and we hope it inspires more research
in this direction.

Figure 2: Overview. We first _render_ an input image to a new camera view
using the disparity. We then _refine_ the image, synthesizing and super-
resolving missing content. As we output both RGB and geometry, this process
can be _repeated_ for perpetual view generation.

##  2 Related Work

Image extrapolation. Our work is inspired by the seminal work of Kaneva et
al.[19], which proposed a non-parametric approach for generating ‘infinite’
images by means of stitching 2D-transformed images, and by patch-based non-
parametric approaches for image extension [29, 1]. We revisit the ‘infinite
images’ concept in a learning framework that also reasons about the 3D
geometry behind each image. Also related to our work are recent deep learning
approaches to the problem of _outpainting_ , i.e., inferring unseen content
outside image boundaries [43, 45, 36], as well as _inpainting_ , the task of
filling in missing content within an image [15, 47]. These approaches use
adversarial frameworks and semantic information for in/outpainting. Our
problem also incorporates aspects of super-resolution [14, 23]. Image-specific
GAN methods also demonstrate a form of image extrapolation and super-
resolution of textures and natural images [50, 34, 30, 33]. In contrast to the
above methods, we reason about the 3D geometry behind each image and study
image extrapolation in the context of temporal image sequence generation.

View synthesis. Many view synthesis methods operate by interpolating between
multiple views of a scene [24, 3, 25, 12, 7], although recent work can
generate new views from just a single input image, as in our work [5, 38, 26,
37, 31, 6]. However, in both settings, most methods only allow for a very
limited range of output viewpoints. Even methods that explicitly allow for
view extrapolation (not just interpolation) typically restrict the camera
motion to small regions around a reference view [49, 35, 8].

One factor that limits camera motion is that many methods construct a static
scene representation, such as a layered depth image [38, 32], multiplane image
[49, 37], or point cloud [26, 44], and inpaint disoccluded regions. Such
representations can allow for fast rendering, but the range of viable new
camera positions is limited by the finite bounds of the scene representation.
Some methods augment this scene representation paradigm, enabling a limited
increase in the range of output views. Niklaus et al. perform inpainting after
rendering (then project back into a point cloud) [26], while SynSin uses a
post-rendering refinement network to produce realistic images from feature
point-clouds [44]. We take inspiration from these methods by rendering and
then refining our output. In contrast, however, our system does not construct
a single 3D representation of a scene. Instead we proceed iteratively,
generating each output view from the previous one, and producing a geometric
scene representation in the form of a disparity map for each frame.

Some methods use video as training data. Monocular depth can be learned from
3D movie left-right camera pairs [22] or from video sequences analysed with
structure-from-motion techniques [4]. Video can also be directly used for view
synthesis [37, 44]. These methods use pairs of images, whereas our model is
trained on sequences of several widely-spaced frames since we want to generate
long-range video.

Video synthesis. Our work is related to methods that generate a video sequence
from one or more images [41, 11, 42, 10, 39, 46]. Many such approaches have
focused on predicting the future of dynamic objects with a static camera,
often using simple videos of humans walking [2] or robot arms [11]. In
contrast, we focus on mostly static scenes with a moving camera, using real
aerial videos of nature. Some recent research addresses video synthesis from
in-the-wild videos with moving cameras [9, 40], but without taking geometry
explicitly into account, and with strict limits on the the length of the
generated video. By accounting for geometry, the task of video prediction
becomes conceptually easier as the movement of pixels from camera motion can
be explicitly modeled using 3D geometry. Our work takes a step towards jointly
modeling the scene geometry and camera motion for video generation.

Figure 3: Illustration of the rendering and refinement steps. (Left): Our
differentiable rendering stage takes a paired RGB image and disparity map from
viewpoint $P_{0}$ and creates a textured mesh representation, which we render
from a new viewpoint $P_{1}$, warping the textures, adjusting disparities, and
returning a binary mask representing regions to fill in. (Right) The
refinement network takes the output of the renderer and uses SPADE [27] as our
network architecture to fill in holes and add details. The output is a new RGB
image and disparity map that can be supervised with reconstruction and
adversarial losses.

##  3 Perpetual View Generation

We introduce _perpetual view generation_ , the task of continually generating
novel views of a scene corresponding to an arbitrary long camera trajectory.
Specifically, at test time, given an RGB image $I_{0}$ and a camera trajectory
$\\{P_{0},P_{1},P_{2},\dots\\}$ of arbitrary length, the task is to output a
new image sequence $\\{I_{0},I_{1},I_{2},\dots\\}$ that forms a video
depicting a flythrough of the scene captured by the initial view. The
trajectory is a series of 3D camera poses
$P_{t}=\left(\begin{smallmatrix}R^{3\times 3}&\mathbf{t}^{3\times 1}\\\
0&1\end{smallmatrix}\right),$ where $R$ and $\mathbf{t}$ are 3D rotations and
translations, respectively. In addition, each camera has an intrinsic matrix
$K$. At test time the camera trajectory may be pre-specified or controlled by
an auto-flight algorithm. At training time camera data is obtained from video
clips via structure-from-motion as in [49].

###  3.1 Approach: Render, Refine, Repeat

We decompose _perpetual view generation_ into three steps, as illustrated in
Fig. 2:

  1. 1.

Render a new view from an old view, by warping the image according to a
disparity map using a differentiable renderer,

  2. 2.

Refine the rendered view and geometry to fill in missing content and add
detail where necessary,

  3. 3.

Repeat this process, generating each future view from the previous one.

Our approach has several desirable characteristics. Representing geometry with
a disparity map allows much of the heavy-lifting of moving pixels from one
frame to the next to be handled by differentiable rendering, ensuring local
temporal consistency. The synthesis task is then reduced to one of image
refinement, which comprises: 1) inpainting disoccluded regions 2) outpainting
of new image regions and 3) super-resolving image content. Here, techniques
from recent image synthesis and generative models can be applied to produce
realistic images. Because every step is fully differentiable, we can train our
refinement network by backpropagating through several view generation
iterations. And because the framework is auto-regressive, novel views may be
infinitely generated with explicit view control, even though training data is
finite in length.

More formally, for an image $I_{t}$ at camera $P_{t}$ we have an associated
disparity map (i.e., an inverse depth map) $D_{t}\in\mathbb{R}^{H\times W}$,
and we compute the next frame $I_{t+1}$ and its disparity $D_{t+1}$ as

$\displaystyle\hat{I}_{t+1},\hat{D}_{t+1},\hat{M}_{t+1}$
$\displaystyle=\mathcal{R}(I_{t},D_{t},P_{t},P_{t+1}),$ (1) $\displaystyle
I_{t+1},D_{t+1}$
$\displaystyle=g_{\theta}(\hat{I}_{t+1},\hat{D}_{t+1},\hat{M}_{t+1}).$ (2)

Here, $\hat{I}_{t+1}$ and $\hat{D}_{t+1}$ are the result of rendering the
image $I_{t}$ and disparity $D_{t}$ from the new camera $P_{t+1}$, using a
differentiable renderer $\mathcal{R}$ [13]. This function also returns a mask
$\hat{M}_{t+1}$ which indicates which regions of the image are missing and
need to be filled in. The refinement network $g_{\theta}$ then inpaints,
outpaints and super-resolves these inputs to produce the next frame $I_{t+1}$
and its disparity $D_{t+1}$. The process is repeated iteratively for $T$ steps
during training, and at test time for an arbitrary length camera trajectory.
Next we discuss each step in detail.

Geometry and Rendering. Our render step $\mathcal{R}$ uses a differentiable
mesh renderer [13]. First, we convert each pixel coordinate $(u,v)$ in $I_{t}$
and its corresponding disparity $d$ in $D_{t}$ into a 3D point in the camera
coordinate system: $(x,y,z)=K^{-1}(u,v,1)/d$. We then convert the image into a
3D triangular mesh where each pixel is treated as a vertex connected to its
neighbors, ready for rendering.

To avoid stretched triangle artefacts at depth discontinuities, and to aid our
refinement network by identifying regions to be completed, we compute a per-
pixel binary mask $M_{t}\in\mathbb{R}^{H\times W}$ by thresholding the
gradient of the disparity image $\nabla\hat{D_{t}}$, computed with a a Sobel
filter:

$M_{t}=\begin{cases}0&\text{where }||\nabla\hat{D_{t}}||>\alpha,\\\
1&\text{otherwise.}\end{cases}$ (3)

The 3D mesh, textured with the image $I_{t}$ and mask $M_{t}$, is then
rendered from the new view $P_{t+1}$, and the rendered image is multiplied
element-wise by the rendered mask to give $\hat{I}_{t+1}$. The renderer also
outputs a depth map as seen from the new camera, which we invert and multiply
by the rendered mask to obtain $\hat{D}_{t+1}$. This use of the mask ensures
that any regions in $\hat{I}_{t+1}$ and $\hat{D}_{t+1}$ that were occluded in
$I_{t}$ are masked out and set to zero (along with regions that were outside
the field of view of the previous camera). These areas are ones that the
refinement step will have to inpaint (or outpaint). See Fig. 2 and Fig. 3 for
examples of missing regions shown in pink.

Refinement and Synthesis. Given the rendered image $\hat{I}_{t+1}$, its
disparity $\hat{D}_{t+1}$ and its mask $M_{t+1}$, our next task is to refine
this image, which encapsulates inpainting, outpainting, and super-resolution.
For this refinement stage we want a model that is both generative and
conditional, so we adopt the state-of-the-art conditional image synthesis
framework of Park et al., which uses spatially-adaptive normalization (SPADE)
to condition the image synthesis process at multiple scales [27]. While the
original SPADE approach generated an image from a semantic segmentation input,
in our work the input is the rendered image, disparity, and mask. The
generator output is a 4-channel image comprising RGB and disparity channels.
We also train a single encoder that encodes the initial input image $I_{0}$ to
compute the latent noise. This module, shown on right in Fig. 3, is trained
with both reconstruction and adversarial losses.

Rinse and Repeat. A crucial part of our approach is to not just refine the RGB
pixels, but also the disparity as well. Together the geometry (represented by
a disparity) and RGB texture provide the necessary information for our
renderer to produce the next view. This insight is key for allowing our
approach to repeat itself perpetually, as each frame and disparity depend only
on the previous one.

Because our render-refinement steps are entirely self-contained, there is no
global scene representation—indeed the only representation is the image and
disparity output from the refinement network at each step. While this is
advantageous because a global representation is expensive to store
perpetually, it means that our render-refine-repeat loop is memory-less and as
a result there is no guarantee of global consistency across multiple
iterations.

Geometric Grounding to Prevent Drift. A notorious challenge in long generation
of sequences is dealing with the accumulation of errors [28]. In a system
where current prediction affects future outputs, subtle deviations in each
iteration may compound, eventually generating predictions that are unseen
during training and causing unexpected behaviors. Repeating the generation
loop in the training process and feeding the network with its own output
ameliorates the drift as can be seen in the _‘No Repeat’_ ablation study
(Section 6). However, we notice that the output, particularly the disparity,
can still drift at test time. Therefore we propose an explicit geometric re-
grounding of the disparity maps.

Specifically, we take advantage of the fact that the rendering process
provides the correct range of disparity from a new viewpoint $\hat{D}_{t+1}$
for visible regions of the previous frame. The refinement network may modify
these values as it refines the holes and blurry regions, which can lead to
drift as the overall disparity gradually may change. However, we can
geometrically correct this by rescaling the refined disparity map to the
correct range by computing a scale factor $\gamma$ via solving

$\min_{\gamma}||M\odot(\log(\gamma D_{t+1})-\log(\hat{D}_{t+1}))||.$ (4)

By scaling the refined disparity by $\gamma$, our approach ensures that the
disparity map stays at a consistent scale, which significantly reduces drift
at test time as shown in Section 6.2.

Figure 4: Processing video for ACID. We run structure from motion on coastline
drone footage collected from YouTube to create the Aerial Coastline Imagery
Dataset (ACID). See Section 4.

##  4 Aerial Coastline Imagery Dataset (ACID)

Learning to generate long sequences requires real image sequences for
training. Many existing datasets for view synthesis do not use sequences, but
only a set of views from slightly different camera positions. Those that do
have sequences are limited in length: RealEstate10K, for example, has
primarily indoor scenes with limited camera movement [49]. To obtain long
sequences with a moving camera and few dynamic objects, we turn to aerial
footage of beautiful nature scenes that are available on the Internet. Nature
scenes are a good starting point for attempting our challenging problem, as
GANs have shown promising results on nature textures [30, 33]. We collected
765 videos using keywords such as ‘coastal’ and ‘aerial footage’, and
processed these videos with SLAM and structure-from-motion following the
approach of Zhou et al. [49], yielding over 13,000 sequences with a total of
2.1 million frames. We make the list of videos and the SfM camera trajectories
available. See Fig. 4 for an illustrative example of our SfM pipeline running
on a coastline video.

Disparity We use the off-the-shelf MiDaS single-view depth prediction method
[22] to obtain disparity maps for every frame. We find that MiDaS is quite
robust and produces sufficiently accurate disparity maps for our method.
Because MiDaS disparity is only predicted up to scale and shift, it must first
be rescaled to match our data. To achieve this, we use the sparse point-cloud
computed for each scene during structure from motion. For each frame we
consider only the points that were tracked in that frame, and apply least-
squares to compute the optimal scale and shift which minimize the disparity
error on these points. We apply this scale and shift to the MiDaS output to
obtain disparity maps $\\{D_{i}\\}$ which are scale-consistent with the SfM
camera trajectories $\\{P_{i}\\}$ for each sequence.

Aligning Camera Speed. The speed of camera motion varies widely in our
collected videos, so we compute a proxy of camera speed in order to normalize
the amount of motion present in training image sequences. We use the
translation magnitude of the estimated camera poses between frames after
scale-normalizing the video as in Zhou et al. [49] to determine a range of
rates at which each sequence can be subsampled in order to obtain a camera
speed within a desired target range. We randomly select frame rates within
this range to subsample videos. We picked a target speed range for training
sequences that varies by up to 30% and, on average, leaves 90% of an image’s
content visible in the next sampled frame. Fig. 4 shows an example of
subsampling.

##  5 Experimental Setup

Losses. We train our approach on a collection of image sequences
$\\{I_{t}\\}_{t=0}^{T}$ with corresponding camera poses
$\\{P_{t}\\}_{t=0}^{T}$ and disparity maps for each frame
$\\{D_{t}\\}_{t=0}^{T}$. Following the literature on conditional generative
models, we use an L1 reconstruction loss on RGB and disparity, a VGG
perceptual loss on RGB [18] and a hinge-based adversarial loss with a
discriminator [27] for the $T$ frames that we synthesize during training. We
also use a KL-divergence loss [21] on our encoder
$\mathcal{L}_{\text{KLD}}=\mathcal{D}_{\text{KL}}(q(z|x)||\mathcal{N}(0,1))$.
Our complete loss function is

$\mathcal{L}=\mathcal{L}_{\text{reconst}}+\mathcal{L}_{\text{perceptual}}+\mathcal{L}_{\text{adversarial}}+\mathcal{L}_{\text{KLD}}$
(5)

The loss is computed over all iterations and over all samples in the mini-
batch.

Method | LPIPS $\downarrow$ | MSE $\downarrow$ | FID $\downarrow$
---|---|---|---
SVG-LP [10] | 0.60 | 0.020 | 135.9
SynSin [44] | 0.32 | 0.018 | 98.1
SynSin-Iter | 0.40 | 0.021 | 143.6
MPI [37] | 0.35 | 0.019 | 65.0
MPI-Iter | 0.47 | 0.020 | 201.2
3D Photos [32] | 0.30 | 0.020 | 123.6
Ours (no-repeat) | 0.30 | 0.022 | 95.4
Ours | 0.32 | 0.020 | 50.6

Table 1: Quantitative evaluation. For LPIPS and MSE we compute the error over
ten frames of ground truth. We report FID-50 computed over all 50 frames
generated from an input test images. See Section 6.1.

Figure 5: FID over time. Left: FID-20 over time for 50 frames generated by
each method. Right: FID-50 over _500_ frames generated by our method using
autopilot. For comparison, we plot FID-50 for the baselines on the first 50
steps. Despite generating sequences an order of magnitude longer, our FID-50
is still lower than that of the baselines. See Sections 6.1, 6.2.

Metrics. Evaluating the quality of the generated images in a way that
correlates with human judgement is a challenge. We use the Fréchet inception
distance (FID), a common metric used in evaluating generative models of
images. FID computes the difference between the mean and covariance of the
embedding of real and fake images through a pretrained Inception network [17]
to measure the realism of the generated images as well as their diversity. We
precompute real statistics using 20k real image samples from our dataset. To
measure changes in generated quality over time, we report FID over a sliding
window: we write FID-$w$ at $t$ to indicate a FID value computed over all
image outputs within a temporal window of width $w$ centered at time $t$, i.e.
$\\{I_{i}\\}$ for $t-w/2<i\leq t+w/2$. For short-range generations where
ground truth images are available, we report mean squared error (MSE) and
LPIPS [48], a perceptual similarity metric that correlates better with human
perceptual judgments than traditional metrics such as PSNR and SSIM.

Implementation Details. We train our model with $T=5$ steps of render-refine-
repeat at an image resolution of 160 $\times$ 256 (as most aerial videos have
a 16:9 aspect ratio). The choice of $T$ is limited by both memory and
available training sequence lengths. The refinement network architecture is
the same as that of SPADE generator in [27], and we also employ the same
multi-scale discriminator. We implement our models in TensorFlow, and train
with a batch size of 4 over 10 GPUs for 7M iterations, which takes about 8
days. We then identify the model checkpoint with the best FID score over a
validation set.

Figure 6: Qualitative comparison over time. We show generated sequence results
for each method at different time steps. Note that we only have ground truth
images for 10 frames; the subsequent frames are generated using an
extrapolated trajectory. Pink region in Ours no-refine indicate missing
content uncovered by the moving camera. Figure 7:  Qualitative comparison. We
show the diversity and quality of many generated results for each method on
the t=15 and 30 frame generation. Competing approaches result in missing or
unrealistic frames. Our approach is able to generate plausible views of the
scene. Figure 8:  Long trajectory generation. From a single image, our
approach can generate 500 frames of video without suffering visually. Please
see the supplementary video for the full effect.

##  6 Evaluation

We compare our approach with three recent state-of-the-art single-image view
synthesis methods—the 3D Photography method of Shih et al. (henceforward ‘3D
Photos’) [32], SynSin [44], and single-view MPIs [37]—as well as the SVG-LP
video synthesis method [10]. We retrain each method on our ACID training data,
with the exception of 3D Photos which is trained on in-the-wild imagery and,
like our method, takes MiDaS disparity as an input. SynSin and single-view MPI
were trained at a resolution of $256\times 256$. SVG-LP takes two input frames
for context, and operates at a lower resolution of $128\times 128$.

The view synthesis methods were not designed for long camera trajectories, so
we also consider iterative variants (SynSin-Iter, MPI-Iter) of these methods
in which instead of synthesizing all output views from the initial input
image, the next viewpoint is produced by using the previously generated output
as the new input at test time. We omit this for 3D photos as iterating is
unfortunately prohibitively slow.

###  6.1 Short-to-medium range view synthesis

To evaluate short-to-medium-range synthesis, we select sequences from our ACID
test data with an input frame and 10 subsequent ground truth frames
(subsampling as described in Section 5), with the camera moving forwards at an
angle of up to 45∘. Although our method is trained on all types of camera
motions, this forward motion is appropriate for comparison with view synthesis
methods which are not designed to handle extreme camera movements.

We then extrapolate the camera motion from the last two frames of each
sequence to extend the trajectory for an additional 40 frames. To avoid the
camera colliding with the scene, we check the final camera position against
the disparity map of the last ground-truth frame, and discard sequences in
which the final pose is outside the image or at a depth large enough to be
occluded by the scene.

This yields a set of 279 sequences with camera trajectories of 50 steps and
ground truth images for the first 10 steps. For short-range evaluation, we
compare to ground truth on the first 10 steps. For medium-range evaluation, we
compute FID scores over all 50 frames.

We apply each method to these sequences to generate novel views corresponding
to the camera poses in each sequence (SVG-LP is the exception in that it does
not take account of camera pose.) See results in Table 5. While our goal is
perpetual view generation, we find that our approach is competitive with
recent view synthesis approaches for short-range synthesis on LPIPS and MSE
metrics. For mid-range evaluation, we report FID-$50$ over all generated 50
frames. Our approach has a dramatically lower FID-50 score than other methods,
reflecting the more naturalistic look of its output. To quantify the
degradation of each method over time, we report a sliding window FID-$20$
computed from $t=10$ to 40. As shown in Fig. 5 (left), the image quality
(measured by FID-$20$) of the baseline methods deteriorates quicker with
increasing $t$ compared to our approach.

Qualitative comparisons of these methods are shown in Fig. 6 and our
supplementary video, which illustrates how the quality of each method’s output
changes over time. Notable here are SVG-LP’s blurriness and inability to
predict any camera motion at all; the increasingly stretched textures of 3D
Photos’ output; and the way the MPI-based method’s individual layers become
noticeable. SynSin does the best job of generating plausible texture, but
still produces holes after a while and does not add new detail. Neither SynSin
nor MPI benefits from being applied iteratively. These results are not
surprising, but illustrate that none of these systems were designed for more
than short-range range synthesis. Fig. 7 shows additional qualitative output
of generating 15 and 30 frames for each method on a variety of inputs.

Ablations. We investigate the benefit of training over multiple iterations of
our render-refine-repeat loop by also training a ‘No Repeat’ variant of our
model with $T=1$. That is, at training time this variant predicts only one
frame ahead instead of the 5 predicted in our full model. At inference time,
the performance on short-range generation, as measured in LPIPS and MSE, is
similar to our full model. But when we look at FID, we observe that this
method generates lower quality images (Table 5), and that they get
substantially worse with increasing $t$ (Fig. 5). This shows the importance of
auto-regressive training to our method.

We next illustrate the contribution of our refine step, via a comparison with
a version of our method which omits this step completely, shown as ‘Ours no-
refine’ in Fig. 6. For clarity, in this figure we set masked pixels to pink at
each step. Note that with increasing $t$, a larger and larger portion of the
image consists of such pixels. In the full model, this region would have been
inpainted or outpainted by our refinement network in one of its preceding
steps. Note also that even non-masked areas of the image are much blurrier
when the refinement step is omitted. This shows the need of our refinement
network in super-resolving image content.

###  6.2 Perpetual view generation

We also evaluate the ability of our model to perform perpetual view generation
by synthesizing videos of 500 frames, using an auto-pilot algorithm to create
an online camera trajectory that avoids flying directly into the ground, sky
or obstacles such as mountains. This algorithm works iteratively in tandem
with image generation to control the camera based on heuristics which measure
the proportion of sky and of foreground obstacles in the scene. See the
supplementary for details. We generate 500 frames for each of our test
sequences and compute their FID-50 over time, as shown in Fig. 5. Our
performance on this metric is robust: even after 500 frames, the FID is lower
than that of all the baseline methods over 50 frames. Fig. 5 also shows the
benefit of our proposed Geometric Grounding—when it’s omitted, the image
quality gradually deteriorates, indicating that drift is an important issue to
resolve.

Fig. 8 shows a qualititative example of long sequence generation. Despite the
challenging problem, our approach retains some aesthetic look of a coastline,
generating new islands, rocks, beaches, and waves as it flies through the
world. The auto-pilot algorithm can receive additional inputs (such as a user-
specified trajectory or random elements), allowing us to generate diverse
videos from a single image. Please see the supplementary video for more
examples and the full effect of these generated fly-through videos.

##  7 Discussion

We introduce a new problem of perpetual view generation and present a novel
framework that combines both geometric and generative techniques to tackle it.
Our system can generate video sequences covering over hundreds of steps, which
to our knowledge has not been shown for prior video or view synthesis methods.
Our results indicate that our hybrid approach is a promising
step—nevertheless, there remain many exciting challenges in this domain.
First, the refinement network, like other GANs, can produce images that seem
realistic but not recognizable [16]. Further advancement in image and video
synthesis generation methods that incorporate geometry would be an interesting
future direction. Second, we have modeled scene generation as a Markov
process. While local temporal consistency is provided by the rendering
process, the memory-less property leads to a system that does not have global
consistency. After some steps, a field may gradually turn into an ocean as we
fly over it. Incorporating memory in this system would be an exciting
direction. Lastly, we do not model dynamic scenes (although our model does a
reasonable job in generating plausible-looking ocean waves): combining our
geometry-aware approach with methods that can reason about object dynamics is
an exciting direction for future work.

## References

  * [1] Connelly Barnes, Eli Shechtman, Adam Finkelstein, and Dan B Goldman.  PatchMatch: A randomized correspondence algorithm for structural image editing.  ACM Transactions on Graphics (Proc. SIGGRAPH), 28(3), Aug. 2009\.
  * [2] Moshe Blank, Lena Gorelick, Eli Shechtman, Michal Irani, and Ronen Basri.  Actions as space-time shapes.  In ICCV, pages 1395–1402. IEEE, 2005.
  * [3] Gaurav Chaurasia, Sylvain Duchêne, Olga Sorkine-Hornung, and George Drettakis.  Depth synthesis and local warps for plausible image-based navigation.  Trans. on Graphics, 32:30:1–30:12, 2013.
  * [4] Weifeng Chen, Shengyi Qian, and Jia Deng.  Learning single-image depth from videos using quality assessment networks.  In The IEEE Conference on Computer Vision and Pattern Recognition (CVPR), June 2019.
  * [5] Xu Chen, Jie Song, and Otmar Hilliges.  Monocular neural image based rendering with continuous view control.  In ICCV, 2019.
  * [6] Xu Chen, Jie Song, and Otmar Hilliges.  Monocular neural image based rendering with continuous view control.  In ICCV, pages 4090–4100, 2019.
  * [7] Inchang Choi, Orazio Gallo, Alejandro Troccoli, Min H Kim, and Jan Kautz.  Extreme view synthesis.  In Proceedings of the IEEE International Conference on Computer Vision, pages 7781–7790, 2019.
  * [8] Inchang Choi, Orazio Gallo, Alejandro Troccoli, Min H Kim, and Jan Kautz.  Extreme view synthesis.  In ICCV, pages 7781–7790, 2019.
  * [9] Aidan Clark, Jeff Donahue, and Karen Simonyan.  Efficient video generation on complex datasets.  arXiv preprint arXiv:1907.06571, 2019.
  * [10] Emily Denton and Rob Fergus.  Stochastic video generation with a learned prior.  arXiv preprint arXiv:1802.07687, 2018.
  * [11] Chelsea Finn, Ian Goodfellow, and Sergey Levine.  Unsupervised learning for physical interaction through video prediction.  In NeurIPS, pages 64–72, 2016.
  * [12] John Flynn, Michael Broxton, Paul Debevec, Matthew DuVall, Graham Fyffe, Ryan Overbeck, Noah Snavely, and Richard Tucker.  Deepview: View synthesis with learned gradient descent.  In The IEEE Conference on Computer Vision and Pattern Recognition (CVPR), June 2019.
  * [13] Kyle Genova, Forrester Cole, Aaron Maschinot, Aaron Sarna, Daniel Vlasic, and William T. Freeman.  Unsupervised training for 3d morphable model regression.  In The IEEE Conference on Computer Vision and Pattern Recognition (CVPR), June 2018.
  * [14] Daniel Glasner, Shai Bagon, and Michal Irani.  Super-resolution from a single image.  In ICCV, pages 349–356, 2009.
  * [15] James Hays and Alexei A Efros.  Scene completion using millions of photographs.  ACM Transactions on Graphics (TOG), 26(3):4–es, 2007.
  * [16] Aaron Hertzmann.  Visual indeterminacy in generative neural art.  arXiv preprint arXiv:1910.04639, 2019.
  * [17] Martin Heusel, Hubert Ramsauer, Thomas Unterthiner, Bernhard Nessler, and Sepp Hochreiter.  Gans trained by a two time-scale update rule converge to a local nash equilibrium.  In NeurIPS, pages 6626–6637, 2017.
  * [18] Justin Johnson, Alexandre Alahi, and Li Fei-Fei.  Perceptual losses for real-time style transfer and super-resolution.  In European conference on computer vision, pages 694–711. Springer, 2016.
  * [19] Biliana Kaneva, Josef Sivic, Antonio Torralba, Shai Avidan, and William T. Freeman.  Infinite images: Creating and exploring a large photorealistic virtual space.  In Proceedings of the IEEE, 2010.
  * [20] Tero Karras, Samuli Laine, and Timo Aila.  A style-based generator architecture for generative adversarial networks.  In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pages 4401–4410, 2019.
  * [21] Diederik P Kingma and Max Welling.  Auto-encoding variational bayes.  arXiv preprint arXiv:1312.6114, 2013.
  * [22] Katrin Lasinger, René Ranftl, Konrad Schindler, and Vladlen Koltun.  Towards robust monocular depth estimation: Mixing datasets for zero-shot cross-dataset transfer.  arXiv preprint arXiv:1907.01341, 2019.
  * [23] Christian Ledig, Lucas Theis, Ferenc Huszár, Jose Caballero, Andrew Cunningham, Alejandro Acosta, Andrew P Aitken, Alykhan Tejani, Johannes Totz, Zehan Wang, et al.  Photo-realistic single image super-resolution using a generative adversarial network.  In CVPR, 2017.
  * [24] Marc Levoy and Pat Hanrahan.  Light field rendering.  In Proceedings of SIGGRAPH 96, Annual Conference Series, 1996.
  * [25] Ben Mildenhall, Pratul P. Srinivasan, Rodrigo Ortiz-Cayon, Nima Khademi Kalantari, Ravi Ramamoorthi, Ren Ng, and Abhishek Kar.  Local light field fusion: Practical view synthesis with prescriptive sampling guidelines.  ACM Transactions on Graphics (TOG), 2019.
  * [26] Simon Niklaus, Long Mai, Jimei Yang, and Feng Liu.  3D Ken Burns effect from a single image.  ACM Transactions on Graphics (TOG), 2019.
  * [27] Taesung Park, Ming-Yu Liu, Ting-Chun Wang, and Jun-Yan Zhu.  Semantic image synthesis with spatially-adaptive normalization.  In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, 2019.
  * [28] Stéphane Ross, Geoffrey Gordon, and Drew Bagnell.  A reduction of imitation learning and structured prediction to no-regret online learning.  In Proceedings of the fourteenth international conference on artificial intelligence and statistics, pages 627–635, 2011.
  * [29] Arno Schödl, Richard Szeliski, David H Salesin, and Irfan Essa.  Video textures.  In Proceedings of the 27th annual conference on Computer graphics and interactive techniques, pages 489–498, 2000.
  * [30] Tamar Rott Shaham, Tali Dekel, and Tomer Michaeli.  Singan: Learning a generative model from a single natural image.  In ICCV, pages 4570–4580, 2019.
  * [31] Lixin Shi, Haitham Hassanieh, Abe Davis, Dina Katabi, and Fredo Durand.  Light field reconstruction using sparsity in the continuous fourier domain.  Trans. on Graphics, 34(1):12:1–12:13, Dec. 2014.
  * [32] Meng-Li Shih, Shih-Yang Su, Johannes Kopf, and Jia-Bin Huang.  3d photography using context-aware layered depth inpainting.  In IEEE Conference on Computer Vision and Pattern Recognition (CVPR), 2020.
  * [33] Assaf Shocher, Shai Bagon, Phillip Isola, and Michal Irani.  Ingan: Capturing and remapping the" dna" of a natural image.  arXiv preprint arXiv:1812.00231, 2018.
  * [34] Assaf Shocher, Nadav Cohen, and Michal Irani.  “zero-shot” super-resolution using deep internal learning.  In CVPR, pages 3118–3126, 2018.
  * [35] Pratul P. Srinivasan, Richard Tucker, Jonathan T. Barron, Ravi Ramamoorthi, Ren Ng, and Noah Snavely.  Pushing the boundaries of view extrapolation with multiplane images.  In The IEEE Conference on Computer Vision and Pattern Recognition (CVPR), June 2019.
  * [36] Piotr Teterwak, Aaron Sarna, Dilip Krishnan, Aaron Maschinot, David Belanger, Ce Liu, and William T Freeman.  Boundless: Generative adversarial networks for image extension.  In Proceedings of the IEEE International Conference on Computer Vision, pages 10521–10530, 2019.
  * [37] Richard Tucker and Noah Snavely.  Single-view view synthesis with multiplane images.  In The IEEE Conference on Computer Vision and Pattern Recognition (CVPR), June 2020.
  * [38] Shubham Tulsiani, Richard Tucker, and Noah Snavely.  Layer-structured 3D scene inference via view synthesis.  In The European Conference on Computer Vision (ECCV), September 2018\.
  * [39] Sergey Tulyakov, Ming-Yu Liu, Xiaodong Yang, and Jan Kautz.  Mocogan: Decomposing motion and content for video generation.  In CVPR, pages 1526–1535, 2018.
  * [40] Ruben Villegas, Arkanath Pathak, Harini Kannan, Dumitru Erhan, Quoc V Le, and Honglak Lee.  High fidelity video prediction with large stochastic recurrent neural networks.  In Advances in Neural Information Processing Systems, pages 81–91, 2019.
  * [41] Carl Vondrick, Hamed Pirsiavash, and Antonio Torralba.  Generating videos with scene dynamics.  In NeurIPS, pages 613–621, 2016.
  * [42] Carl Vondrick and Antonio Torralba.  Generating the future with adversarial transformers.  In CVPR, pages 1020–1028, 2017.
  * [43] Yi Wang, Xin Tao, Xiaoyong Shen, and Jiaya Jia.  Wide-context semantic image extrapolation.  In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pages 1399–1408, 2019.
  * [44] Olivia Wiles, Georgia Gkioxari, Richard Szeliski, and Justin Johnson.  SynSin: End-to-end view synthesis from a single image.  In CVPR, 2020.
  * [45] Zongxin Yang, Jian Dong, Ping Liu, Yi Yang, and Shuicheng Yan.  Very long natural scenery image prediction by outpainting.  In Proceedings of the IEEE International Conference on Computer Vision, pages 10561–10570, 2019.
  * [46] Yufei Ye, Maneesh Singh, Abhinav Gupta, and Shubham Tulsiani.  Compositional video prediction.  In ICCV, 2019.
  * [47] Jiahui Yu, Zhe Lin, Jimei Yang, Xiaohui Shen, Xin Lu, and Thomas S Huang.  Free-form image inpainting with gated convolution.  In ICCV, 2019.
  * [48] Richard Zhang, Phillip Isola, Alexei A. Efros, Eli Shechtman, and Oliver Wang.  The unreasonable effectiveness of deep features as a perceptual metric.  In CVPR, 2018.
  * [49] Tinghui Zhou, Richard Tucker, John Flynn, Graham Fyffe, and Noah Snavely.  Stereo magnification: Learning view synthesis using multiplane images.  ACM Trans. Graph., 37(4):65:1–65:12, 2018.
  * [50] Yang Zhou, Zhen Zhu, Xiang Bai, Dani Lischinski, Daniel Cohen-Or, and Hui Huang.  Non-stationary texture synthesis by adversarial expansion.  arXiv preprint arXiv:1805.04487, 2018.


Principal Component Analysis and Hidden Markov Model for Forecasting Stock
Returns

Eugene W. Park

A thesis submitted in partial fulfillment
of the requirements for the degree of
Master of Science
Courant Institute of Mathematical Sciences
New York University
May, 2023

###### Abstract

This paper presents a method for predicting stock returns using principal
component analysis (PCA) and the hidden Markov model (HMM), and tests the
results of trading stocks based on this approach. Principal component analysis
is applied to the covariance matrix of stock returns for companies listed in
the S&P 500 index, and interpreting principal components as factor returns, we
apply the HMM model on them. Then we use the transition probability matrix and
state conditional means to forecast the factors returns. Reverting the factor
returns forecasts to stock returns using eigenvectors, we obtain forecasts for
the stock returns. We find that, with the right hyperparameters, our model
yields a strategy that outperforms the buy-and-hold strategy in terms of the
annualized Sharpe ratio.

Keywords: Principal component analysis, factor model, hidden Markov model,
stock market, forecasting

##  1 Introduction

Stock market forecasting has been a prolonged practice of interest for people
in various fields of discipline, and thus, various approaches have been used
to tackle the problem from classical time series analysis and using factor
models to high frequency trading and using deep learning techniques.

A widely used method to analyze time series data, the hidden Markov model has
been a popular method to analyze financial markets. [1], [2], [4], and [5]
have used the HMM on its own to make predictions of stock prices, while others
have combined the HMM with other methods such as the long short term memory
model [3] and fuzzy logic [6]. The usage of the hidden Markov model on stock
price prediction is more extensive than what has been cited, but within the
author’s knowledge, there seems to be more attention on using HMM as a
preliminary step than combining techniques that refine the input for the HMM.

In this paper, we present a method that applies principal component analysis
(PCA) as a form of a factor model to preprocess our data and uses the HMM on
the preprocessed data to forecast of stock returns. Then, we test the accuracy
of this model by computing the directional accuracy of the forecasts, and
evaluate the long-term, risk-adjusted returns of various trading strategies
based on this model.

The rest of this paper is organized as follows. In section 2, we give a brief
overview of the PCA as a factor model, and in section 3, briefly introduce the
HMM. Then, we explain our model in detail in section 4, and discuss the
implementation and results of our model as trading strategies in section 5. We
conclude the paper in section 6 by mentioning some shortcomings of our paper
and providing ways for improvement.

##  2 Brief Introduction: Principal Component Analysis as Factor Models

###  2.1 Principal Component Analysis

Given a dataset $X\in\mathbb{R}^{t\times n}$ where the rows of $X$ are de-
meaned samples of the random vector $\tilde{X}\in\mathbb{R}^{n}$, principal
component analysis is performed as follows:

  * $\diamond$

Compute the covariance matrix of $X$:

$M=X^{T}X\in\mathbb{R}^{n\times n}$ (1)

  * $\diamond$

Compute the eigendecomposition of $M$, which exists by the spectral theorem
for symmetric matrices, such that

$M=EGE^{T}$

where $E\in\mathbb{R}^{n\times n}$ is the column matrix of n orthonormal
eigenvectors $u_{1},\dots,u_{n}\in\mathbb{R}^{n}$, and
$G\in\mathbb{R}^{n\times n}$ is a diagonal matrix of eigenvalues
$\lambda_{1},\dots,\lambda_{n}$.

Since we are allowed to reorder the eigenvalues in descending order with the
eigenvectors ordered accordingly, we may assume that

$\lambda_{1}>\dots>\lambda_{n}$

Each $\lambda_{i}$ denotes the amount of variance explained in $M$, and
$u_{i}$ are the corresponding directions called the ”principal directions.”

Note that $u_{1}$ is the direction of the highest variance in $M$ and $u_{n}$
is the direction of the lowest variance.

The ”principal components” are then the dataset $X$ in the directions of the
principal direction $u_{i}$:

$w_{i}=\begin{bmatrix}w_{i,1}\\\ \vdots\\\ w_{i,n}\end{bmatrix}=Xu_{i}$

###  2.2 Factor Models

The factor model, introduced in the Arbitrage Pricing Theory (APT) in [8],
states that there exist explanatory variables (called ”factors”) that explain
the systematic behavior of asset returns. Hence,

$\textbf{r}=B\textbf{f}+\epsilon$ (2)

where

  * $\diamond$

$\textbf{r}\in\mathbb{R}^{n\times 1}$ is a random vector of returns for n
assets,

  * $\diamond$

$\textbf{f}\in\mathbb{R}^{k\times 1}$ is a random vector of returns for the k
factors with $\mathbb{E}(\textbf{f})=\mu_{f}$ and $\sigma^{2}(\textbf{f})=F$,

  * $\diamond$

$B\in\mathbb{R}^{n\times k}$ consists of columns representing factor loadings
for the k factors,

  * $\diamond$

$\epsilon\in\mathbb{R}^{n\times 1}$ is random noise for which we assume
$\mathbb{E}(\epsilon)=0$, $\sigma^{2}(\epsilon)=D$, a diagonal matrix.

From , 2 we get

$\mathbb{E}(\textbf{r})=B\mu_{\textbf{f}}$
$\sigma^{2}(\textbf{r})=\Sigma_{\textbf{r}}=BFB^{T}+D$ (3)

where
$\Sigma_{\textbf{r}}=\mathbb{E}(\textbf{r}^{T}\textbf{r})\in\mathbb{R}^{n\times
n}$ is the covariance matrix of r.

###  2.3 Classical Factor Models and PCA Factor Models

***This subsection is largely based on [7] Suppose we have a dataset $X$ of n
asset returns across t time ($X\in\mathbb{R}^{t\times n}$). The classical,
APT-type, factor models uses features of the assets to derive the k factors
that explain $X$. Hence, factors are exogenous to $X$. The Fama-French five
factor model [9] and factor-based risk models built by MSCI are some of the
popular examples of the classic factor models.

The PCA factor models differ from the classical ones in that the factors are
not discerned exogenously. The factors are data-driven (i.e. inferred from the
data matrix $X$) by applying the eigendecompoistion on the covariance matrix
in equation 1:

$M=X^{T}X=EGE^{T}$

The PCA factor model sets the eigenvectors as factor loadings,

$E=B$ (4)

, and the factor returns as the n principal components

$\textbf{f}=XE=[w_{1},\dots w_{n}]$ (5)

Then, we have

$F=\textbf{f}^{T}\textbf{f}=(XE)^{T}(XE)$ $\Rightarrow
EFE^{T}=E(XE)^{T}(XE)E^{T}=X^{T}X$ $\Rightarrow F=G$

Assuming the PCA model to be exact, hence leaving aside $D$ in (3) (i.e.
$\epsilon=0$ in (2)),

$M=X^{T}X=EGE^{T}=BFB^{T}+0$ (6)

Thus, the PCA factor model defines n factors that explain the full covariance
matrix of our dataset of returns, $X$, where each factor corresponds to $X$
rotated in n independent direction $u_{i}\in E$, $i=1,\dots,n$, and the factor
loadings are these n independent directions. †††Note that the independence of
eigenvectors come from the spectral theorem of symmetric matrices, and that
the independence of $u_{i}$’s suggests that $f|_{k}$ consists of independent
column vectors The analogue to (2) is then

$X^{T}=E(XE)^{T}$ (7)

where $X^{T}\in\mathbb{R}^{n\times t}$ and $(XE)^{T}\in\mathbb{R}^{n\times t}$
consists of $t$ samples in the columns.

We note that the although factors derived by PCA are less intuitive, the PCA
offers the benefit of discovering factors that not may be discovered
exogenously.

##  3 Brief Introduction: Hidden Markov Model

The fundamental argument of the hidden Markov model (HMM) is that underlying
the observed sequence of time series data $Y=\\{y_{1},\dots,y_{t}\\}$, there
exists a Markov chain, $Z=\\{z_{1},\dots,z_{t}\\}$ that generates $Y$ where
each $z_{i}\in S=\\{1,...,N\\}$, the state space. The underlying Markov chain
has an initial distribution $\textbf{L}=\\{L_{1},\dots,L_{n}\\}$ and a
transition probability matrix $\textbf{P}=(p_{ij})_{i,j\in S}$, where $p_{ij}$
denotes the probability of $z_{l}$ transitioning from state i to j. Finally,
assuming our observations to be in state space $O$, the emission probability
matrix, $\textbf{R}=(r_{ij})_{i\in S,j\in O}$, denotes the probability of
observation given that we are in a certain state:
$r_{ij}=\mathbb{P}(Y=j|Z=i)$.

Hence, the parameters that define the HMM are:

$\Theta=(\textbf{L},\textbf{P},\textbf{R})$ (8)

and we are mainly concerned with the following problems:

  1. 1.

Selecting the best model $\Theta$ given a range of model options, $\Theta_{k}$
and the sequence of observations $Y$:

$\underset{\Theta_{k}}{\text{argmax}}\mathbb{P}(Y|\Theta_{k})$

  2. 2.

Determining the most probable state sequence $Z$ given $Y$ and $\Theta$:

$\underset{Z}{\text{argmax}}\mathbb{P}(Z|Y,\Theta)$

  3. 3.

Estimating the parameters $\Theta$ given $Y$:

$\underset{\Theta}{\text{argmax}}\mathbb{P}(Y|\Theta)$

The forward algorithm is an algorithm to solve the first problem, Viterbi
algorithm for the second, and Baum-Welch Algorithm or EM (expectation-
maximization) algorithm for the last.‡‡‡Refer to [10] for details of the
algorithms

##  4 PCA + HMM as a Forecasting Model

In this section, we discuss details for how we use PCA with HMM to develop a
model that forecasts asset returns.

Suppose that we are working with a dataset, $X$, of returns for n assets
through T time periods such that the columns of $X$ are time series data of
return for a particular asset:

$X=(r_{ij})_{i=\\{1,\dots,T\\},j=\\{1,\dots,n\\}}$

, where $r_{ij}$ is the return of company j at time i. Thus, our goal is to
forecast the returns on the n assets for the next period,
$\hat{X}_{T+1,n=\\{1,\dots,n\\}}$. Furthermore, since the directional accuracy
of our forecasts is crucial for our trading strategies, we are, in fact,
primarily concerned with:

$\text{sign}(\hat{X}_{T+1,n=\\{1,\dots,n\\}})$ (9)

###  4.1 Implementing PCA

PCA only requires our dataset to be de-meaned, but we normalize $X$ for each
column and denote it as $Y$:

$Y=\frac{X-\mu_{x}}{\sigma_{X}}$ (10)

We apply PCA by computing the eigendecomposition of the covariance matrix of
$Y$ such that the eigenvalues are in decreasing order and the eigenvectors are
ordered accordingly:

$H=Y^{T}Y=EGE^{T},\ H,E,G\in\mathbb{R}^{n\times n}$
$\text{diag}G=\\{\lambda_{1},\dots,\lambda_{n}\\},\
\lambda_{1}>\dots>\lambda_{n}$

Then, as in (7), we have a full PCA factor model:

$Y^{T}=E(YE)^{T}$

We assume that the covariance matrix $H$, or equivalently $Y$, contains some
noise. Thus, we want to de-noise our dataset before training the HMM.

Recall that the eigenvalues suggest the percentage of variance (i.e.
information) of $H$ explained by rotating it in the direction of the
corresponding eigenvectors. Suppose the amount of noise in the covariance
matrix is $p\%$. Then, we take the first $k$ eigenvectors that explain at
least $(1-p)\%$ of variance such that

$\lambda_{1}+\dots+\lambda_{k}\geq 1-p$ $\lambda_{1}+\dots+\lambda_{k-1}<1-p$

and compute the corresponding principal components, ($w_{1},\dots,w_{k}$),
which represent k of the n factor returns as in (5).

Then, we extract $\approx p\%$ §§§$\approx p\%$ is due to the fact that the k
eigenvectors does not necessarily, in fact, most often does not, explain
exactly $(1-p\%)$ of the matrix. of noise from $H$ by restricting our set of
eigenvectors, $E$, to the first k eigenvectors

$E|_{k}=\\{u_{1},\dots,u_{k}\\}\in\mathbb{R}^{n\times k}$

and computing k factor returns,

$f|_{k}=\\{w_{1},\dots,w_{k}\\}=YE|_{k}=\begin{bmatrix}y_{1}\cdot
u_{1}&\dots&y_{1}\cdot u_{k}\\\ y_{2}\cdot u_{1}&\dots&y_{2}\cdot u_{k}\\\
\vdots&\vdots&\vdots\\\ y_{T}\cdot u_{1}&\dots&y_{T}\cdot
u_{k}\end{bmatrix}\in\mathbb{R}^{T\times k}$ (11)

, where $y_{i}\in\mathbb{R}^{k}$ is the returns for the k factors in time $i$.
Note that $f|_{k}$ consists of columns that represent $Y$ in k independent
directions with the most amount of information (i.e. $f|_{k}$ consists of
independent column vectors).

Assuming the full PCA factor model to contain some noise, we are essentially
assuming that not all of the principal component factors defined by PCA are
significant factors; factors that explain small portions of variance of the
covaraince matrix are just noise. Hence, we have:

$Y^{T}=E(YE)^{T}=E|_{k}(YE|_{k})^{T}+\epsilon_{k}$ (12)

where $\epsilon_{k}$ is the error term that contains $\approx p\%$ of
information in $H$.

###  4.2 HMM on Factor Returns

Now, we extract the noise from $Y$ and apply the HMM to $f|_{k}$, the time
series of $k$ factor returns, to obtain a one-step-ahead forecast for these
factor returns. Then, we revert the forecast of the factors back to assets
which will then be our forecast for the assets.

Before, we begin training the model, we assume the emission probability in
each state in the state space, $S$, to have have a Gaussian distribution.
Hence, the parameters of our model are

$\Theta=(\textbf{L},\textbf{P},\textbf{R})$
$\Rightarrow\Theta=(\textbf{L},\textbf{P},\mu_{i},\sigma^{2}_{i}),i=1,\dots,N$
(13)

, where N is the number of states in the state space $S$, $\mu_{i}$ is the
mean of the observation sequence in state i, and $\sigma^{2}_{i}$ is the
variance of observation sequence in state i.

####  4.2.1 Training & Model Selection

To determine the number of state space, $N$, we train the Gaussian HMM (13)
for each $j$ in [2,3,4,5,6,7,8], denoted as $\Theta_{j}$, compute the log
likelihood, $\ln(\mathbb{P}(Y|\Theta_{j}))$, using the forward algorithm, and
compute Akaike information criterion (AIC)[11] for each model.

$\text{AIC}=-2\ln(\ln(\mathbb{P}(Y|\Theta_{j})))+2j$

Then, we choose our model to be $\Theta_{j}$, which yields the lowest AIC.

In the training step of each $\Theta_{j}$, we use the Baum-Welch algorithm to
estimate the parameters, and thus, choose the initial conditions for our
parameters as follows:

$\textbf{L}_{0}=1/j$ (14) $\textbf{P}_{0}=1/j^{2}$ (15)
$\mu_{i,0}=\mathbb{E}(o)_{i}$ (16) $\sigma^{2}_{i,0}=\text{var}(o)_{i}$ (17)

, where $o$ is the observation sequence, $o=\\{o_{1},\dots,o_{T}\\}$. Here, we
are assuming that it is equally likely for the underlying Markov chain to
start in any state (14) and that it is equally likely for the Markov chain to
transit from one state to another (15). (16) suggests that we assume the
$\mu_{i}$ to be the sample mean of our observations for all states
$i=1,\dots,j$ (i.e. same initial mean for all states), and (17) suggests that
we assume $\sigma^{2}_{i}$ to be the sample variance of our observations for
all states $i=1,\dots,j$ (i.e. same initial variance for all states).

Note that since $f|_{k}$ consists of k time series of factor returns, our
observation sequnce is k-dimensional.

$o=\\{o_{1},\dots,o_{T}\\},\text{where }o_{t}=\\{o_{t}^{1},\dots,o_{t}^{k}\\}$

Hence, $\sigma^{2}_{i}$ is a $k\times k$ covariance matrix. Since, the k
factor returns are independent, however, we may assume the covariance matrix
to be diagonal. Thus, in (17), $\sigma^{2}_{i,0}=\text{var}(o)_{i}$ is a
$k\times k$ diagonal covariance matrix, for state $i$, with sample variances
of $o=\\{o^{1},\dots,o^{k}\\}$ in the diagonal.

####  4.2.2 Forecasting

Now, once we have found and trained our model, which we will denote again as
$\Theta$, we use the estimated parameters to forecast the returns of the k
factors in the next time period, $(\hat{f|_{k}})_{T+1}$, as follows:

  1. 1.

We use the Viterbi algorithm to obtain the sequence of the underlying Markov
chain:

$Z=\\{z_{1},\dots,z_{T}\\},\ z_{i}\in S=\\{1,\dots,N\\}$

, where $z_{T}=i$ states that the underlying Markov chain is currently in
state $i$.

  2. 2.

Then, we use the transition probability matrix, $\hat{\textbf{P}}$, and the
state conditional mean, $\hat{\mu}_{i}$, to obtain an estimate for
$(f|_{k})_{T+1}$:

$(\hat{f}|_{k})_{T+1}=\displaystyle\sum_{j\in S}^{N}P_{ij}\hat{\mu}_{j}$ (18)

Since $\hat{\mu}_{j}\in\mathbb{R}^{k}$, for all $j=\\{1,\dots,N\\}$,
$(\hat{f}|_{k})_{T+1}\in\mathbb{R}^{k}$.

Thus, our forecast of the k factor returns in the next period depends on 1)
how likely it is to transition from state $i$ to $j$, and 2) estimated state
conditional means of the k factors.

Although we may trade factors, especially if it were defined exogenously as in
the classical factor models, it is hard to trade endogenous factors defined by
the PCA. Hence, we revert the forecasts for the factor returns to forecasts
for the $n$ assets. I.e. we want:

$(\hat{f}|_{k})_{T+1}\rightarrow\hat{Y}_{T+1}$

Recall from (12) that

$Y^{T}=E(YE)^{T}=E(f)^{T}$

where $f=\\{w_{1},\dots,w_{n}\\}$. Since we only have an estimate of $f|_{k}$
for time $T+1$, we first estimate $f_{k+1},\dots,f_{n}$ for $T+1$. Recall that
the best estimate is its mean. Thus, we have, for $j=k+1,\dots,n$,

$(\hat{f}_{j})_{T+1}=\mathbb{E}((f_{j})_{T+1})=\frac{1}{T}\displaystyle\sum_{i=0}^{T}(Y_{i}\cdot
u_{j})=\frac{1}{T}\displaystyle\sum_{i=0}^{T}(Y_{i})\cdot u_{j}=0$ (19)

where $Y_{i}\in\mathbb{R}^{n}$ is the i-th row of $Y$, and the last equality
follows from $Y$ being normalized with row mean equal to 0. Equivalently,

$\mathbb{E}((\hat{\epsilon_{k}})_{T+1})=\mathbb{E}(\hat{Y}_{T+1}^{T}-E|_{k}(\hat{Y}_{T+1}E|_{k})^{T})=\frac{1}{T}\displaystyle\sum_{h=0}^{T}Y_{h}^{T}-E|_{k}(\frac{1}{T}\displaystyle\sum_{h=0}^{T}Y_{h}^{T})E|_{k})^{T}=0$

where $Y_{h}^{T}$ are the h-th column of $Y^{T}$.

Now, we have

$\hat{f}_{T+1}=\\{(\hat{f}|_{k})_{T+1},0,\dots,0\\}$

and from (12), we have

$\hat{Y}^{T}_{T+1}=E(\hat{f}_{T+1})^{T}=E|_{k}(\hat{f}|_{k})^{T}_{T+1}$
$\Rightarrow\hat{Y}_{T+1}=(\hat{f}|_{k})_{T+1}E|_{k}^{T}$

where $\hat{Y}_{T+1}\in\mathbb{R}^{1\times n}$,
$(\hat{f}|_{k})_{T+1}\in\mathbb{R}^{1\times k}$, and
$E|_{k}^{T}\in\mathbb{R}^{k\times n}$. Thus, we obtain forecasts for the $n$
assets by multiplying the transpose of the k eigenvectors to the forecasts for
the k factor returns.

Finally, since $Y$ is the normalized version of $X$ (10),

$\hat{X}_{T+1}=\sigma_{X}(\hat{y}_{T+1})+\mu_{X}$ (20)

##  5 Implementation

We use data for the weekly returns of companies listed in the S&P500 to train
our model.¶¶¶Data was obtained from finance.yahoo.com∥∥∥List obtained from
https://en.wikipedia.org/wiki/List_of_S%26P_500_companies on April 28, 2023
The returns are calculated using the closing stock prices, and each return
sequence comprise the columns of our dataset, denoted $X$. Hence, we have
$X\in\mathbb{R}^{Txn}$, where $T$ is the last period of the time span of our
data and $n$ is the number of assets.

We note that, as companies are newly listed and delisted from the S&P500
index, some companies that are listed in the index at the time of data
retrieval may not have been listed previously, and some may have been delisted
during the time span of our data. Thus, we are subject to survivor bias, but
minimize this effect by ensuring that we have at least 400 data points of
company returns for each time period. With at least 80% of the index, we
conjecture that the several companies we fail to capture only explains a
minimal amount of information/variance in the covariance matrix, $X^{T}X$,
which we implement the PCA on. In other words, we suspect that most of the
information from the companies we fail to capture have already been captured
by the eigenvectors corresponding to large eigenvalues, and that the
additional information from those companies constructs eigenvectors with very
small eigenvalues. Thus, we suspect that they are eventually extracted during
the process of denoising our covariance matrix via PCA. Because this is only a
conjecture, however, we do not know for sure how much we are affected by
survivor bias, and thus, leaving room for improvement upon this paper.

###  5.1 Model Training

Ensuring our dataset to include returns data for at least 400 companies listed
in the S&P500 index for each $t$, we start our data from July 23, 2004. Using
10 years’ worth of weekly data, we train our PCA+HMM model and forecast the
weekly returns for the next period. Then, rolling the 10 years window by one
week******The exact time samples (rows of our dataset) for each window may
slightly differ due to different numbers of holidays and trading days every
year. Moreover, due to the entering and exiting of companies from the index
since the beginning of our dataset, the number of stocks (columns of our
dataset) may also differ., we forecast returns for a total of 100 weeks (the
week for which we make our last forecast is July 8, 2016, thus our data spans
from July 23, 2004 to July 8,2016)††††††Additional data is available, but we
limit ourselves to 100 forecasts for computational convenience.. Note that,
since we train the PCA+HMM for each time window, the total number of PCA
factors, the total number of states in the HMM, and the model parameters may
vary over time.

Since PCA is essentially an eigendecomposition, We use the linalg module in
the numpy library in Python to implement PCA. To train the HMM model, we use
the gaussianhmm module in the hmmlearn Python library.

Recall from part 4.1 that we have to specify the hyperparameter, $p$, that
specifies the minimum percentage of noise to be extracted during the dimension
reduction process. We train our model with four choices of $p$:

$p=\\{45\%,30\%,15\%,10\%\\}$

which corresponds to keeping no less than 55%, 70%, 85%, and 90% of the
information in $Y$, our normalzied return matrix.

####  5.1.1 Model Implications

Before, we test trading strategies based on our model, we make several
observations of our trained model.

Figure 1 shows the total number of PCA factors kept in the model throughout
time for the different choices of $p$. Notice that as we keep more noise in
our model, the number of PCA factors kept not only increases, but also changes
more frequently. Moreover, there is also a more gradual change in the total
number of PCA factors in the market, such as from the first 20 periods to the
period between the 20th and 60th week in figures 1.a and 1.b. There is also a
gradual shift in 1.c and 1.d from the first 50 weeks to the next. This
suggests that there may exist short-term PCA factors and longer-term PCA
factors.

(a) P = 45%

(b) P = 30%

(c) P = 15%

(d) P = 10%

Figure 1: Number of PCA Factors

Now, we observe the state-conditional means and variances calibrated from the
HMM model. Recall that the HMM calibrates the means and variances for $K$ PCA
factors for each state, and since we run the model for 100 periods, we have a
series of state-conditional means and variances as follows:

$\mu=\\{\mu_{1},\dots,\mu_{100}\\}$
$\sigma^{2}=\\{\sigma_{1}^{2},\dots,\sigma_{100}^{2}\\}$

where $\mu_{t},\sigma_{t}^{2}\in\mathbb{R}^{N_{t}\times k_{t}}$. Note that
$N_{t}$ and $k_{t}$ has a time subscript since the total number of states and
PCA factors vary.

At each time, $t$, we average $\mu_{t},\sigma_{t}^{2}$ along the k factors to
get a time series of state-conditional means and variances (of returns)

$\tilde{\mu}=\\{\tilde{\mu}_{1},\dots,\tilde{\mu}_{100}\\}$
$\tilde{\sigma}^{2}=\\{\tilde{\sigma}_{1}^{2},\dots,\tilde{\sigma}_{100}^{2}\\}$

where $\tilde{\mu}_{t},\tilde{\sigma}_{t}^{2}\in\mathbb{R}^{N_{t}}$. Moreover,
since the total number of states are changing, we sort $\tilde{\mu}_{t}$ in
decreasing order with $\tilde{\sigma}_{t}^{2}$ having the corresponding order.
Thus, we disregard the actual state numbers, such as rather we are in state 1
or state 2, and put meaning on the states by ranking them according to the
calibrated state-conditional means averaged across the $k_{t}$ PCA factors.
For each $p$, we find the minimum number states that existence throughout time
and visualize the emission probability distributions in figure 2.

(a) P = 45%

(b) P = 30%

(c) P = 15%

(d) P = 10%

Figure 2: State-Conditional Emission Probability Distributions

Note:x-axis in percentages. 3 refers to 300%.

From figure 2, we can see that more noise being kept in our factor returns
yields probability distributions that are less spread across states,
suggesting that the underlying states lose significance and the HMM model
becomes less useful. Furthermore, figure 2.a, the model with $45\%$ of noise
extracted from the normalized return dataset,$Y$, is intuitively
interpretable. State 4 to seems to be ”bear market” state, while state 0 to
can be viewed as the ”bull market” state, and states 1,2, and 3 are the
relatively stable states of the market. We can also view state 4 as where
investors have excessively fearful, state 0 as the state where they are
excessively optimistic, states 1 and 3 as states where they are moderately
optimistic and fearful, and state 2 as where they neutral about the market.
Thus, models with more spread-out probability distributions provide more
reliable and interpretable information about the market which may be useful
for task such as risk-management.

###  5.2 Trading Strategies

Now we construct and test simple trading strategies based on our model.

  1. 1.

Our first trading strategy is solely based on the forecasts given by
LABEL:asset_forecasts. We short the assets for which $\hat{X}_{T+1}$ predicts
a negative return and long the ones for which $\hat{X}_{T+1}$ predicts a
positive return. We weight the $n$ stocks equally (i.e. 1/n) and trade at the
closing time of the stock market at the last trading day of every week.

  2. 2.

The second strategy depends on the forecasts of the normalized asset returns,
$\hat{y}_{T+1}$. Since

$Y=\frac{X-\mu_{X}}{\sigma_{X}}$

, we hypothesize that there is predictive trading signal in not the return
itself, but the excess return over the average.

When testing our trading strategies, we assume the following:

  1. 1.

we are able to trade at exactly the closing time

  2. 2.

we can exactly match the $n$ assets with equal weight in our portfolio

  3. 3.

and we do not incur any trading costs.

###  5.3 Results

To evaluate the performance of our trading strategies, we compare our trading
strategies with the simple buy-and-hold strategy of all companies in the
index. As measurements of performance, we compute:

  1. 1.

the ”winning probability” of our forecasts (i.e. probability of correctly
forecasting the direction of the returns in the next period).

  2. 2.

the annualized Sharpe ratio:

$S=\frac{r-rf}{\sigma}$ (21)

where $r$ is the asset return, $rf$ is the return of the risk-free rate, and
$\sigma$ is the standard deviation of asset. We will assume the risk-free rate
to be zero.

Table 1: Winning Probability

p% | Strategy 1 | Strategy 2
---|---|---
45% | 0.532 | 0.490
30% | 0.543 | 0.504
15% | 0.538 | 0.511
5% | 0.532 | 0.490

  * $\diamond$

Note: Rounded to the nearest thousandth.

The winning probabilities of the two trading strategies are shown in table 1.
We observe that strategy 1 clearly outperforms strategy 2 in forecasting the
direction of the returns. Yet, strategy 1 yield just above 50%, suggesting
that significant profit can be made only in a casino-style type of trading
where the model trades constantly without going bankrupt. Hence, for this
strategy to be profitable over the long-run, additional features must be
included such as the trade size.

Table 2 shows the Sharpe ratios of the trading strategies, and we observe that
the strategy with the highest Sharpe ratio is strategy 2 with $p=15\%$. This
outperforms the buy-&-hold strategy by more than 50%. From figure 5, we
observe that the strategy yields lower cumulative return than the buy-&-hold,
but is more stable; it rides out two ocassions of downturns, during late-2014
and early-2016, without much losses, contrary to the other strategies. It also
has a winning probability of 51%, and thus, seems to be useful in terms of
risk management and as a stable component within a portfolio.

We also note that the Sharpre ratios depend heavily on the hyperparameter $p$.
The difference between the highest and lowest ratios of strategy is as high as
0.91. For strategy 1, there seems to be a trend where the ratio falls as $p$
decreases to $30\%$, but increases again as more noise is extracted. This
suggests the significance of finding the right hyperparameter $p$, which could
yield significantly higher Sharpe ratios than the ones presented here.

Table 2: Annualized Sharpe Ratio

p% | Strategy 1 | Strategy 2
---|---|---
45% | 0.688 | 0.45
30% | 0.581 | 0.581
15% | 0.703 | 1.36
10% | 0.877 | 0.726
Buy-&-Hold | 0.828 |

  * $\diamond$

Note: Values are percentages.

Figure 3: Return Sequences for Strategies 1 and 2 with p = 45% Figure 4:
Return Sequences for Strategies 1 and 2 with p = 30% Figure 5: Return
Sequences for Strategies 1 and 2 with p = 15% Figure 6: Return Sequences for
Strategies 1 and 2 with p = 10%

##  6 Conclusion

We have shown that with the right hyperparameter $p$, the PCA + HMM model
provides an intuitive understanding of the market through the estimated
emission probability distributions of the states, and yields model-based
trading strategies with forecast accuracy slightly above 50%. We have also
seen that these strategies have the potential to outperform the market by
finding the right $p$. This work suggests that future research into finding
the optimal p will be very useful. Moreover, we have made several assumptions
which such as zero transaction costs and that we are able to trade at exactly
the closing prices at the closing time. These assumptions are unrealistic and
fails in practice. Hence, our paper could be improved by evaluating the
performance under realistic conditions, such as existence of slippage,
transactions costs, and possibly price impact as well (when we assume large
trade size). Acknowledgement: We thank professors Kenneth Winston, Petter
Kolm, and Jonathan Goodman for their valuable comments and suggestions
throughout the progress of this work.

## Appendix

Appendix A.1. Proof of Forecast Equation 18:
Denote $r_{t}$, return at time t, as the rows of $f_{t}$.

$\displaystyle(\hat{f}|_{k})_{T+1}$
$\displaystyle=\mathbb{E}(r_{T+1}|Z_{T}=i)$
$\displaystyle=\displaystyle\sum_{j\in
S}^{N}\mathbb{E}(r_{T+1},Z_{T+1}=j|Z_{T}=i)$
$\displaystyle=\displaystyle\sum_{j\in
S}^{N}\mathbb{E}(r_{T+1}|Z_{T+1}=j)\mathbb{P}(Z_{T+1}=j|Z_{T}=i)$
$\displaystyle=\displaystyle\sum_{j\in S}^{N}P_{ij}\hat{\mu}_{j}$

## References

  * [1] N. Nyugen. ”Hidden Markov Model for Stock Trading”. International Journal of Financial Studies, pages 6-36, 2018.
  * [2] B. Dhingra, A. Gupta. ”Stock Market Prediction Using Hidden Markov Models”. 2012 Students Conference on Engineering and Systems, pages 1-4, 2012.
  * [3] J. Huo, M. Liu, Y. Wu, J. Wu. ”Stock Market Trend Analysis Using Hidden Markov Model and Long Short Term Memory”. eprint arXiv:2104.09700, 2012\.
  * [4] R. Hassan, B. Nath. ”Stock market forecasting using hidden Markov model: a new approach”. 5th International Conference on Intelligent Systems Design and Applications, pages 192-196, 2005.
  * [5] G. Kavitha, A. Udhayakumar, D. Nagarajan. ”Stock Market Trend Analysis Using Hidden Markov Models”. International Journal of Computer Science and Information Security , 2013\.
  * [6] R. Hassan. ”A Combination of hidden Markov model and fuzzy model for stock market forecasting”. Neurocomputing, pages 3439-3446, 2009.
  * [7] K. J. Winston. Quantitative Risk and Portfolio Management: Theory and Practice. Cambridge University Press, 2023 (forthcoming).
  * [8] S. A. Ross. ”The arbitrage theory of capital asset pricing”. Journal of Economic Theory,  pages 341-360, 1976.
  * [9] E. F. Fama, K.R. French. ”A five-factor asset pricing model”. Journal of Financial Economics,  pages 1-22, 2014.
  * [10] T. Li, W. E, E. Vanden-Eijnde. Applied Stochastic Analysis. American Mathematical Society, 2019.
  * [11] A. Hirotugu. ”A new look at the statistical model identification”. IEEE TRansactions on Automatic Control, pages 716-723, 1974.