Structural Periodic Vector Autoregressions

Daniel Dzikowski***TU Dortmund University, Department of Statistics, D-44221
Dortmund, Germany<EMAIL_ADDRESS>corresponding author
Carsten Jentsch†††TU Dortmund University, Department of Statistics, D-44221
Dortmund, Germany<EMAIL_ADDRESS>
TU Dortmund University TU Dortmund University

August 28, 2024

###### Abstract

While seasonality inherent to raw macroeconomic data is commonly removed by
seasonal adjustment techniques before it is used for structural inference,
this approach might distort valuable information contained in the data. As an
alternative method to commonly used structural vector autoregressions (SVAR)
for seasonally adjusted macroeconomic data, this paper offers an approach in
which the periodicity of not seasonally adjusted raw data is modeled directly
by structural periodic vector autoregressions (SPVAR) that are based on
periodic vector autoregressions (PVAR) as the reduced form model. In
comparison to a VAR, the PVAR does allow not only for periodically time-
varying intercepts, but also for periodic autoregressive parameters and
innovations variances, respectively. As this larger flexibility leads also to
an increased number of parameters, we propose linearly constrained estimation
techniques. Overall, SPVARs allow to capture seasonal effects and enable a
direct and more refined analysis of seasonal patterns in macroeconomic data,
which can provide useful insights into their dynamics. Moreover, based on such
SPVARs, we propose a general concept for structural impulse response analyses
that takes seasonal patterns directly into account. We provide asymptotic
theory for estimators of periodic reduced form parameters and structural
impulse responses under flexible linear restrictions. Further, for the
construction of confidence intervals, we propose residual-based (seasonal)
bootstrap methods that allow for general forms of seasonalities in the data
and prove its bootstrap consistency. A real data application on industrial
production, inflation and federal funds rate is presented, showing that useful
information about the data structure can be lost when using common seasonal
adjustment methods.

Keywords: Periodic Vector Autoregressions, Structural VARs, Linear
Restrictions,
Seasonal Adjustment, Residual-Based Seasonal Bootstrap

##  1 Introduction

Macroeconomic time series often exhibit seasonal structure defined as
recurrent intra-year movement. The reasons for seasonal fluctuations in
economic time series are manifold. Classical literature as, for instance,
Hylleberg (1992), Hylleberg et al. (1993), Canova & Hansen (1995) and Franses
(1996) argue that seasonality of time series may be caused by calendar and
weather effects or by seasonally varying behavior of economic agents, among
others. The standard approach to deal with seasonality is to remove the
seasonal structure of macroeconomic time series by using seasonal adjustment
methods prior to a further analysis. The most popular approaches are the
classical X-11 seasonal adjustment program established by Shiskin et al.
(1967), the TRAMO/SEATS program by Gómez & Maravall Herrero (1996) and model-
based upgrades X-12-ARIMA (Findley et al., 1998) and X-13-ARIMA-SEATS
(Monsell, 2007). In the literature, seasonal adjustment methods are often
controversial due to their high complexity, lack of linkage to economic theory
and the resulting lack of transparency; see, for instance, Gersovitz &
MacKinnon (1978), Bell & Hillmer (1984), Osborn & Smith (1989) and Franses
(1996). Ghysels et al. (1996) and Ghysels & Osborn (2001) figure out that
seasonal adjustment methods may distort the structure of data leading to
misspecified standard errors and inference. Recently, Doppelt (2021) argue
that the results of SVAR analyses are highly dependent on the used seasonal
adjustment method and structural analyses of not seasonally adjusted data
should be prefered. Furthermore, even after applying seasonal adjustment
methods, the seasonally adjusted data may still display seasonal structure;
see Del Barrio Castro & Osborn (2004), Bell (2012) and Findley et al. (2017).
As an example, Rudebusch et al. (2015) and Lunsford (2017) found residual or
leftower seasonality in seasonally adjusted GDP estimates by the Bureau of
Economic Analysis (BEA). Even after BEA introduced a new estimation strategy
to remove residual seasonality of seasonally adjusted GDP and GDI and their
components in 2018, Consolvo & Lunsford (2019) detected remaining residual
seasonality in adjusted GDP estimates after BEAs improvements.

Direct modeling of not seasonally adjusted macroeconomic data by periodic
autoregressions can be an alternative to VARs fitted after seasonal
adjustment. Periodic autoregressions are useful to capture not only seasonal
structure present in the mean of time series data, but also in their
covariances in the sense that both the mean and the (auto)covariances are
periodic functions of time. Processes with periodically varying covariances,
often called periodically correlated stochastic processes, have a long
history. First introduced by Gladyshev (1961), pioneering work on (univariate)
periodic autoregressive (PAR) processes was done by Jones & Brelsford (1967),
Pagano (1978) and Troutman (1979), who derived inference techniques on
periodic autoregressive parameter estimates, predictions using PAR models and
the connection between PAR and the related stationary multivariate
autoregressive process. PARs for modeling seasonally not adjusted
macroeconomic data have been used by Osborn (1988), Osborn & Smith (1989),
Novales & de Fruto (1997), Franses & Van Dijk (2005) and Ghysels et al.
(2006), who mainly focused on forecasting studies. However, although PAR
models have received some attention in economics and certain methodological
developments have been made (see e.g. Franses & Paap, 2004), they have not
received much popularity for structural analysis. This is mainly because, for
analyzing not seasonally adjusted macroeconomic data, a vector-valued
extension of the PAR model class is required that naturally comes with a lot
of parameters.

Asymptotic theory for such periodic vector autoregressive (PVAR) models have
been developed by Ula (1990, 1993), Franses & Paap (2004) and Aknouche (2007),
who derived periodic stationarity conditions and parameter estimation methods.
Further, P(V)AR models are often considered in stochastic volatility modeling
of financial data, see e.g. Bollerslev & Ghysels (1996), Aknouche & Bibi
(2009) and Aknouche (2017). Thanks to their large number of parameters, PVAR
models are extremely flexible. However, unconstrained PVAR estimation does
often not work well in applications with moderate sample sizes as typically
encountered in macro econometrics due to pronounced overfitting issues. To
address these, Ursu & Duchesne (2009) and Boubacar Maïnassara & Ursu (2023)
derived least squares estimators under linear parameter restrictions. However,
they only allow to impose linear restriction on parameters _within_ a given
season to control the number of parameters to be estimated. In this paper, we
extend their approach by explicitly allowing also for practically linear
constraints _across_ seasons.

Figure 1: Estimated SDs (top panels) and ACFs (bottom panels) of monthly
seasonally not adjusted industrial production (IP), inflation (INF) and
federal funds rate (FFR) from 1968-2019.

The importance of constraints across seasons in a PVAR setup is illustrated in
our real data application, in which we analyze monthly data consisting of not
seasonally adjusted industrial production (IP), CPI inflation (INF) and the
federal funds rate (FFR). In Figure 1, we present estimated spectral densities
(SDs) and autocovariance functions (ACFs) for all three (univariate) time
series. Taking a closer look at the raw data, we find that IP and INF exhibit
strong periodic annual patterns, while FFR shows no periodic structure at all.

In Figure 2, the estimated SDs and ACFs of _seasonally demeaned_ IP, INF and
FFR are displayed. Apparently, while seasonal demeaning removes almost all
periodicity of INF, it is clearly not able to adequately remove the entire
periodic structure of IP. That is, even after seasonally demeaning, IP still
exhibits pronounced periodic structure dominated by sub-year periods. Further,
also in view of the large number of parameters of an unconstraint PVAR, since
two out of three variables have hardly any periodic structure, it is not
necessary to fit a full PVAR model, i.e., without any constraints across
seasons, to explain the dependence structure of the data.

Figure 2: Estimated SDs (top panels) and ACFs (bottom panels) of monthly
seasonally demeaned industrial production (IP), inflation (INF) and federal
funds rate (FFR) from 1968-2019.

This paper is constructed as follows. Section 2 introduces PVAR processes,
their stationary higher-dimensional VAR representation and states conditions
for (periodic) stationarity. In Section 3, we find a representation of PVAR
models suitable for estimation purposes and derive (constrained) least-squares
estimators for the PVAR parameters. In comparison to Ursu & Duchesne (2009)
and Boubacar Maïnassara & Ursu (2023), we extend their approach for the
constrained estimation of PVAR parameters by allowing also for practically
relevant constraints _across seasons_. We provide asymptotic theory for the
(constrained) estimators of the PVAR parameters under (possibly weak) periodic
white noise assumptions. Section 4 introduces the structural PVAR model, which
enables a direct structural analysis of periodic time series without prior
seasonal adjustment, and states asymptotic properties of (structural) impulse
responses. Section 5 proposes suitable residual-based seasonal bootstrap
methods for the construction of confidence intervals and provides bootstrap
consistency results. Following the analysis of Peersman (2005), a real data
application on IP, INF and FFR is provided in Section 6, while Section 7
summarizes our findings. All proofs and additional figures are deferred to
supplementary material.

##  2 Periodic Vector Autoregressive Process

In this section, we discuss periodic vector autoregressive (PVAR) processes
and derive their higher-dimensional time-invariant VAR and moving average
representations. Let $\\{y_{t}\\}_{t\in\mathbb{Z}}$ be an $m$-dimensional PVAR
process (Ursu & Duchesne, 2009) described by the model equation

$\displaystyle
y_{Sn+s}=\nu(s)+A_{1}(s)y_{Sn+s-1}+\dots+A_{p(s)}(s)y_{Sn+s-p(s)}+\epsilon_{Sn+s},\quad
n\in\mathbb{Z},$ (1)

where $S$ is the number of seasons within one cycle and the $y_{Sn+s}$ denotes
the observation in season $s$ and cycle $n$. The model order $p(s)$, the
$m$-dimensional seasonal intercept $\nu(s)$, and the ($m\times m$) seasonal
autoregressive coefficient matrices $A_{i}(s)$, $i=1,\dots,p(s)$ are allowed
to vary across the seasons $s=1,\dots,S$. The innovation process
$\\{\epsilon_{t}\\}_{t\in\mathbb{Z}}$ is assumed to be an $m$-dimensional
_periodic white noise process_ with period $S$ and non-singular covariance
matrices $Var(\epsilon_{Sn+s})=\Sigma_{\epsilon}(s)$, $s=1\dots,S$ meaning
that $\\{\epsilon_{t}\\}_{t\in\mathbb{Z}}$ can be represented as

$\displaystyle\epsilon_{Sn+s}=H_{0}(s)w_{Sn+s},\quad n\in\mathbb{Z},$ (2)

where $\\{w_{t}\\}_{t\in\mathbb{Z}}$ is white noise with
$E(w_{t}w_{t}^{\prime})=\bm{I}_{m}$ and the non-singular matrix $H_{0}(s)$
satisfies $\Sigma_{\epsilon}(s)=H_{0}(s)H_{0}(s)^{\prime}$ for all
$s=1,\dots,S$. Here, $\bm{I}_{m}$ denotes the $m$-dimensional identity matrix
and $\Sigma_{\epsilon}(s)$ is allowed to vary across the seasons
$s=1,\dots,S$. Since the model order $p(s)$ of a PVAR process can vary across
seasons as well, we introduce the notation PVAR($\bm{p}$) with
$\bm{p}=(p(1),\dots,p(S))$ and PVAR($p$), if $p(s)=p$ for all $s=1,\dots,S$.
For $S=1$ the PVAR($\bm{p}$)=PVAR($p$) process collapses to a common VAR($p$)
process.

However, the PVAR($\bm{p}$) in (1) can also be represented with constant order
$p$, where $p=\max\\{p(1),\dots,p(S)\\}$ by imposing zero restrictions on
$A_{j}(s)$ for $p(s)<j\leq p$ (Lund & Basawa, 2000). A higher-dimensional
time-invariant, i.e. stationary, VAR($P$) representation of the PVAR($\bm{p}$)
process in (1) is given by the model equation (Franses & Paap, 2004)

$\displaystyle\bm{A}_{0}Y_{n}=\nu+\bm{A}_{1}Y_{n-1}+\dots+\bm{A}_{P}Y_{n-P}+\xi_{n},\quad
n\in\mathbb{Z},$ (3)

where

$\displaystyle Y_{n}$
$\displaystyle=(y_{Sn+1}^{\prime},y_{Sn+2}^{\prime},\dots,y_{Sn+S}^{\prime})^{\prime},\quad\xi_{n}$
$\displaystyle=(\epsilon_{Sn+1}^{\prime},\epsilon_{Sn+2}^{\prime},\dots,\epsilon_{Sn+S}^{\prime})^{\prime},\quad\nu$
$\displaystyle=(\nu(1)^{\prime},\dots,\nu(S)^{\prime})^{\prime}$

are of dimension $Sm$, respectively. The coefficient matrices
$\bm{A}_{0},\bm{A}_{1}\dots,\bm{A}_{P}$ contain the PVAR coefficient matrices
from (1) and are of dimension ($Sm\times Sm$). Precisely, it holds

$\displaystyle\bm{A}_{0}$
$\displaystyle=\begin{pmatrix}\bm{I}_{m}&\bm{0}&\dots&\bm{0}\\\
\scalebox{0.6}[1.0]{$-$}A_{1}(2)&\bm{I}_{m}&\ddots&\vdots\\\
\vdots&\ddots&\ddots&\bm{0}\\\
\scalebox{0.6}[1.0]{$-$}A_{S\scalebox{0.6}[0.6]{$-$}1}(S)&\dots&\scalebox{0.6}[1.0]{$-$}A_{1}(S)&\bm{I}_{m}\end{pmatrix},\bm{A}_{i}=\begin{pmatrix}A_{Si}(1)&A_{Si\scalebox{0.6}[0.6]{$-$}1}(1)&\dots&A_{Si\scalebox{0.6}[0.6]{$-$}S\scalebox{0.6}[0.6]{$+$}1}(1)\\\
A_{Si\scalebox{0.6}[0.6]{$+$}1}(2)&A_{Si}(2)&\ddots&\vdots\\\
\vdots&\ddots&\ddots&A_{Si\scalebox{0.6}[0.6]{$-$}1}(S\scalebox{0.6}[1.0]{$-$}1)\\\
A_{Si\scalebox{0.6}[0.6]{$+$}S\scalebox{0.6}[0.6]{$-$}1}(S)&\dots&A_{Si\scalebox{0.6}[0.6]{$+$}1}(S)&A_{Si}(S)\end{pmatrix}$

for $i=1,\ldots,P$. The VAR model order $P$ of the process
$\\{Y_{n}\\}_{n\in\mathbb{Z}}$ is given by $P=\lceil p/S\rceil$, where
$\lceil.\rceil$ is the ceiling function. Since the number of unit roots of a
PVAR($\bm{p}$) process equals the one in its higher-dimensional VAR($P$)
representation (Franses & Paap, 2004), the periodic stationarity and causality
condition for a PVAR($\bm{p}$) process is satisfied if and only if its higher-
dimensional VAR($P$) form is stationary and causal. As $\det(\bm{A}_{0})=1$,
according to Lütkepohl (2005), model (3) is equivalent to

$\displaystyle
Y_{n}=\bm{A}_{0}^{-1}\nu+\bm{A}_{0}^{-1}\bm{A}_{1}Y_{n-1}+\dots+\bm{A}_{0}^{-1}\bm{A}_{P}Y_{n-P}+\bm{A}_{0}^{-1}\xi_{n},\quad
n\in\mathbb{Z},$

such that $\\{Y_{n}\\}_{n\in\mathbb{Z}}$ is stationary and causal if and only
if

$\displaystyle\det(\bm{I}_{Sm}-\bm{A}_{0}^{-1}\bm{A}_{1}z-\dots-\bm{A}_{0}^{-1}\bm{A}_{P}z^{P})\neq
0$ (4)

for all $z\in\mathbb{C}$ with $|z|\leq 1$. When the higher-dimensional VAR
representation of a PVAR($\bm{p}$) satisfies the stationarity and causality
condition, $\\{y_{t}\\}_{t\in\mathbb{Z}}$ is (called) periodically stationary
and a causal infinite order moving-average representation of (1) exists, which
is given by

$\displaystyle
y_{Sn+s}=\mu(s)+\sum\limits_{k=0}^{\infty}\Phi_{k}(s)\epsilon_{Sn+s-k},\quad
n\in\mathbb{Z},$ (5)

where $\Phi_{k}(s)$ are $(m\times m)$-dimensional matrices consisting of the
moving average coefficients and $\mu(s)=E(y_{Sn+s})$, $s=1,\dots,S$. If the
PVAR($\bm{p}$) is periodically stationary, its moving average coefficient
matrices $\Phi_{k}(s)$ are absolutely summable for all $s=1,\dots,S$.
Following Ursu & Duchesne (2009), they can be obtained recursively by the
relations

$\displaystyle\Phi_{0}(s)=\bm{I}_{m}\quad\text{and}\quad\Phi_{k}(s)=\sum\limits_{j=1}^{k}A_{j}(s)\Phi_{k-j}(s-j),\quad
k\in\mathbb{N}.$

The moving average coefficient matrices are periodic in the sense that
$\Phi_{k}(Sn+s)=\Phi_{k}(s),\;n\in\mathbb{Z},\;k\in\mathbb{N}_{0}.$ Impulse
responses are defined as responses of the system of variables $k$ time steps
after a one-time unit shock at time $t$ (Kilian & Lütkepohl, 2017), while
moving average coefficient matrices can be interpreted as the effect of a one-
time unit shock at time $t-k$ on $y_{t}$. For a standard VAR, the impulse
responses are given by its moving average coefficient matrices due to the
time-invariant structure of VARs. However, in PVARs, this equality does not
hold since the shocks’ effect on the system of variables depends on the season
the shock occurs. Hence, in order to get PVAR impulse responses
$\Phi_{k}^{IR}(s)$, we need to shift the seasonal index of periodic moving
average coefficient matrices $k$-steps further, that is

$\Phi_{k}^{IR}(s)=\Phi_{k}(s+k)=\sum\limits_{j=1}^{k}A_{j}(s+k)\Phi_{k-j}(s+k-j),\;\;\;k\in\mathbb{N}_{0}.$

For $i,j=1,\dots,m$, the element in the $i$-th row and $j$-th column of
$\Phi_{k}^{IR}(s)$ denotes the impulse response of the $i$-th component of
$y_{Sn+s+k}$ to a unit shock in the $j$-th component of $\epsilon_{Sn+s}$.
Consequently, the PVAR impulse responses are also periodic through the direct
link to the periodic moving average coefficient matrices. There is also a
direct link between PVAR impulse responses and the impulse responses of its
higher dimensional VAR representation. If the stationarity condition (4)
holds, the VAR($P$) in (3) has an MA representation

$\displaystyle
Y_{n}=\bm{A}_{0}^{-1}\nu+\sum\limits_{h=0}^{\infty}\Pi_{h}\bm{A}_{0}^{-1}\xi_{n-h},\quad
n\in\mathbb{Z},$

where the moving average coefficient matrices $\Pi_{h}$ of the VAR process
$\\{Y_{n}\\}_{n\in\mathbb{Z}}$ defined by (3) can be represented as
$\Pi_{h}=J\bm{A}^{h}J^{\prime}$, $h\in\mathbb{N}_{0}$ with $(Sm\times SmP)$
matrix $J=[\bm{I}_{Sm},\bm{0},\dots,\bm{0}]$ and the VAR($P$) companion matrix
$\bm{A}$ is of dimension $(SmP\times SmP)$. The VAR($P$) companion matrix
$\bm{A}$ is defined as the parameter matrix of the higher-dimensional VAR(1)
representation of the VAR($P$). The VAR impulse responses $\Pi_{h}^{IR}$ are
then obtained by

$\displaystyle\Pi_{h}^{IR}=\Pi_{h}\bm{A}_{0}^{-1}=\begin{pmatrix}\Phi_{Sh}^{IR}(1)&\Phi_{Sh-1}^{IR}(2)&\dots&\Phi_{Sh-S+1}^{IR}(S)\\\
\Phi_{Sh+1}^{IR}(1)&\Phi_{Sh}^{IR}(2)&\ddots&\vdots\\\
\vdots&\ddots&\ddots&\Phi_{Sh-1}^{IR}(S)\\\
\Phi_{Sh+S-1}^{IR}(1)&\dots&\Phi_{Sh+1}^{IR}(S-1)&\Phi_{Sh}^{IR}(S)\end{pmatrix},\;\;h\in\mathbb{N}_{0}.$
(6)

##  3 Estimation of PVARs and Asymptotic Inference

In this section, we provide asymptotic theory for a least squares estimation
approach for the PVAR model parameters, which allows for flexible linear
restrictions. In this regard, we extend the setup of Ursu & Duchesne (2009)
and cover more general forms of linear restrictions, which also allow for
constraints across seasons111In comparison to Ursu & Duchesne (2009), we find
a representation of the PVAR($\bm{p}$) process that allows for imposing
general linear restrictions on the parameter vector. Ursu & Duchesne (2009)
consider each season of the PVAR($\bm{p}$) process individually and are
therefore only able to restrict the seasonal parameters separately such that
no constraints across seasons can be imposed. In our approach, this is
possible and enables, for instance, to constrain a PVAR($p$) to a VAR($p$) for
$S$ greater than one.. One may think, it is easier to use the time-invariant
VAR representation of a PVAR to derive asymptotic properties of the PVAR
estimators, but the VAR representation needs to be inverted which makes it
more complicated than deriving the asymptotic results directly using the
individual PVAR components (Ursu & Duchesne, 2009).

In addition, it should be noted that unconstrained PVAR estimation does not
work well in most applications due to pronounced overfitting issues.
Consequently, general linear restrictions play a crucial role in our setup. To
address these issues, we derived a PVAR setup which is also able to collapse
to a non-periodic VAR under suitable linear restrictions.

###  3.1 Estimation of PVAR Models

In order to estimate the PVAR parameters suitably under general linear
restrictions, a multivariate least squares approach is used. We assume that we
observe $y_{1},\dots,y_{T}$ and, for notational convenience, suppose pre-
sample values $y_{s^{\prime}},\dots,y_{0}$, where
$s^{\prime}=\min\\{1-p(1),\dots,S-p(S)\\}$ are available and the sample size
$T=SN$ can be divided into $N$ (complete) cycles of $S$ seasons per cycle.
Then, the PVAR($\bm{p}$) process in (1) can be rewritten as

$\displaystyle Z=BX+E,$ (7)

where $Z,E,B$ and $X$ are defined as

$\displaystyle Z$
$\displaystyle=(y_{1},\dots,y_{SN}),\;\;\;\;\;\;E=(\epsilon_{1},\dots,\epsilon_{SN}),$
$\displaystyle B$
$\displaystyle=\bigl{(}\nu(1),A_{1}(1),\dots,A_{p(1)}(1),\nu(2),A_{1}(2),\dots,A_{p(2)}(2),\dots,\nu(S),A_{1}(S),\dots,A_{p(S)}(S)\bigr{)},$
$\displaystyle X$
$\displaystyle=\begin{pmatrix}X_{0}(1)&&0&X_{1}(1)&&0&&X_{N-1}(1)&&0\\\
&\ddots&&&\ddots&&\dots&&\ddots&\\\
0&&X_{0}(S)&0&&X_{1}(S)&&0&&X_{N-1}(S)\end{pmatrix},$

which are of dimension $(m\times SN)$, $(m\times SN)$,
$(m\times\sum_{s=1}^{S}(mp(s)+1))$ and $(\sum_{s=1}^{S}(mp(s)+1)\times SN)$,
respectively. The periodic autoregressive parameters are included in $B$,
while the entries $X_{n}(s)$ of $X$ are $(mp(s)+1)$-dimensional random vectors
given by
$X_{n}(s)=(1,y_{Sn+s-1}^{\prime},\dots,y_{Sn+s-p(s)}^{\prime})^{\prime}$.
Vectorizing the model equation (7) yields

$\displaystyle z=vec\\{Z\\}$
$\displaystyle=vec\\{BX\\}+vec\\{E\\}=\\{X^{\prime}\otimes I_{m}\\}\beta+e,$
(8)

where $z=(y_{1}^{\prime},\dots,y_{SN}^{\prime})^{\prime}$, $\beta=vec\\{B\\}$
and
$e=vec\\{E\\}=(\epsilon_{1}^{\prime},\dots,\epsilon_{SN}^{\prime})^{\prime}$,
respectively. Note that $\beta$ is given by
$\beta=(\beta^{\prime}(1),\dots,\beta^{\prime}(S))^{\prime}$ with
$\beta(s)=vec\\{(\nu(s),A_{1}(s),\dots,A_{p(s)}(s))\\}$. Since the covariance
matrix of $\epsilon_{Sn+s}$ is given by $\Sigma_{\epsilon}(s)$ and
$\\{w_{t}\\}_{t\in\mathbb{Z}}$ in (2) is an (uncorrelated) white noise
process, the covariance matrix $\Sigma_{E}$ of $e$ is block-diagonal with
$\Sigma_{E}=I_{N}\otimes\Sigma_{\epsilon}$, where
$\Sigma_{\epsilon}=\text{diag}\bigl{(}\Sigma_{\epsilon}(1),\Sigma_{\epsilon}(2),\dots,\Sigma_{\epsilon}(S)\bigr{)}\in\mathbb{R}^{Sm\times
Sm}$ is also block-diagonal.

Ursu & Duchesne (2009) and Boubacar Maïnassara & Ursu (2023) use linear
restrictions imposed on parameters _within_ one season and do not allow for
restrictions _across_ seasons. However, restrictions across seasons can play a
crucial role, especially in macroeconomic applications, since they allow
parameters to be constrained as constant across the seasons. In general, any
set of linear constraints imposed on $\beta$ can be represented by

$\displaystyle\beta=R\gamma+r,$ (9)

where $R$ is a known $(m\sum_{s=1}^{S}(mp(s)+1)\times M)$-dimensional matrix
with rank $M$, $r$ is a $m\sum_{s=1}^{S}(mp(s)+1)$-dimensional known vector
and $\gamma$ is an $M$-dimensional unconstrained parameter vector of interest.
Substituting $\beta=R\gamma+r$ into (8) yields

$\displaystyle\mathbb{z}_{r}=\\{X^{\prime}\otimes\bm{I}_{m}\\}R\gamma+e,$

where $\mathbb{z}_{r}=z-\\{X^{\prime}\otimes\bm{I}_{m}\\}r$. Hence, the least
squares estimator $\widehat{\gamma}_{LS}$ of $\gamma$ is defined as the
minimizer of the residual sum of squares $RSS(\gamma)$, where

$\displaystyle RSS(\gamma)=\sum\limits_{t=1}^{SN}\
\epsilon_{t}^{\prime}\epsilon_{t}=e^{\prime}e=[\mathbb{z}_{r}-\\{X^{\prime}\otimes\bm{I}_{m})\\}R\gamma]^{\prime}[\mathbb{z}_{r}-\\{X^{\prime}\otimes\bm{I}_{m})\\}R\gamma].$

Hence, according to Lütkepohl (2005), the corresponding least squares
estimator is given by
$\widehat{\gamma}=\left[R^{\prime}\\{XX^{\prime}\otimes\bm{I}_{m}\\}R\right]^{-1}R^{\prime}\\{X\otimes\bm{I}_{m}\\}\mathbb{z}_{r}$
and the restricted least squares estimator $\widehat{\beta}_{res}$ for $\beta$
is obtained by replacing $\gamma$ by $\widehat{\gamma}$ in $\beta=R\gamma+r$,
that is

$\displaystyle\widehat{\beta}_{res}=R\widehat{\gamma}+r=R\left[R^{\prime}\\{XX^{\prime}\otimes\bm{I}_{m}\\}R\right]^{-1}R^{\prime}\\{X\otimes\bm{I}_{m}\\}\mathbb{z}_{r}+r.$
(10)

The unrestricted least-squares estimator $\widehat{\beta}_{LS}$ is achieved by
setting $R=\bm{I}_{m\sum_{s=1}^{S}(mp(s)+1)}$ and $r=\bm{0}$ in (10) leading
to
$\widehat{\beta}_{LS}=\left[\\{XX^{\prime}\otimes\bm{I}_{m}\\}\right]^{-1}\\{X\otimes\bm{I}_{m}\\}z$.
On the other hand, if $p(s)=p$ for all $s=1,\dots,S$ and $R$ and $r$ are
defined as $R=\bm{1}_{S}\otimes\bm{I}_{m(mp+1)}$ and $r=\bm{0}$, the seasonal
PVAR($p$) estimators $\beta(s)$ are restricted to be the same across seasons,
where $\bm{1}_{l}$ denotes the $l$-dimensional vector of ones. Once the
periodic autoregressive parameters are estimated for all $s=1,\dots,S$, a
natural candidate to estimate the periodic covariance matrices
$\Sigma_{\epsilon}(s)$ is

$\widehat{\Sigma}_{\epsilon}(s)=\frac{1}{N-k(s)}\sum\limits_{n=0}^{N-1}\widehat{\epsilon}_{Sn+s}\widehat{\epsilon}_{Sn+s}^{\prime},$

where $k(s)$ denotes the number freely varying parameters in season $s$. In
the following, let $\sigma(s)=vech\\{\Sigma_{\epsilon}(s)\\}$ and denote by
$\widehat{\sigma}(s)=vech\\{\widehat{\Sigma}_{\epsilon}(s)\\}$ its empirical
counterpart, where the $vech$ operator stacks the entries on the lower-
triangular part of a square matrix columnwise below each other. Further, let
$\sigma=(\sigma(1)^{\prime},\dots,\sigma(S)^{\prime})^{\prime}$ and
$\widehat{\sigma}=(\widehat{\sigma}(1)^{\prime},\dots,\widehat{\sigma}(S)^{\prime})^{\prime}$,
which are both of dimension $S\widetilde{m}$ with $\widetilde{m}=m(m+1)/2$.

###  3.2 Asymptotic Theory for PVAR Estimators

In order to derive the joint asymptotic properties of $\widehat{\beta}_{res}$
under general linear restrictions and periodic covariance matrix estimators
$\widehat{\sigma}$, some notation is introduced. At first, for $s=1\dots,S$
and $k\in\mathbb{N}$, we define $mp(s)\times m$ matrices
$C_{k}(s)=(\Phi_{k-1}^{\prime}(s-1),\ldots,\Phi_{k-p(s)}^{\prime}(s-p(s))^{\prime}$,
where $\Phi_{k}(s)$ denote the moving average coefficient matrices in (5). In
addition, let $L_{mS}$ be the duplication matrix that satisfies
$(vech\\{B_{1}\\}^{\prime},\dots,vech\\{B_{S}\\}^{\prime})^{\prime}=L_{mS}(vec\\{B_{1}\\}^{\prime},\dots,vec\\{B_{S}\\}^{\prime})^{\prime}$
for symmetric ($m\times m$) matrices $B_{1},\dots,B_{S}$. Further, for
$a,b,c\in\mathbb{Z}$ and $s_{1},s_{2}=1,\ldots,S$, we define the ($m^{2}\times
m)$ and ($m^{2}\times m^{2}$) dimensional matrices

$\displaystyle\kappa_{a,b}(s_{1},s_{2})\;=E(vec\\{\epsilon_{Sn+s_{1}}\epsilon_{Sn+s_{1}-a}^{\prime}\\}\epsilon_{S(n-b)+s_{2}}^{\prime}),$
(11)
$\displaystyle\tau_{a,b,c}(s_{1},s_{2})=\begin{cases}E(vec\\{\epsilon_{Sn+s_{1}}\epsilon_{Sn+s_{1}-a}^{\prime}\\}vec\\{\epsilon_{S(n-b)+s_{2}}\epsilon_{S(n-b)+s_{2}-c}^{\prime}\\})&\\\
\quad\quad-
vec\\{\Sigma_{\epsilon}(s_{1})\\}vec\\{\Sigma_{\epsilon}(s_{2})\\}^{\prime},&a=c=0\\\
E(vec\\{\epsilon_{Sn+s_{1}}\epsilon_{Sn+s_{1}-a}^{\prime}\\}vec\\{\epsilon_{S(n-b)+s_{2}}\epsilon_{S(n-b)+s_{2}-c}^{\prime}\\}),&\text{
else}\end{cases}$ (12)

and the $(Sm^{2}\times Sm^{2})$ matrices
$\tau_{a,b,c}=(\tau_{a,b,c}(s_{1},s_{2}))_{s_{1},s_{2}=1,\dots,S}$. The
symbols $\overset{d}{\to}$ and $\overset{p}{\to}$ denote convergence in
distribution and probability, respectively, and
$\mathcal{N}(\bm{\mu},\bm{\Sigma})$ denotes a normal distribution with mean
vector $\bm{\mu}$ and covariance matrix $\bm{\Sigma}$.

For deriving asymptotic theory for weak PVARs, we impose the following
assumptions on the process $\\{y_{t}\\}_{t\in\mathbb{Z}}$, which allows the
PVAR innovations to be a periodic weak white noise.

######  Assumption 3.1 (PVAR with Periodic Weak White Noise).

  * (i)

Let the true DGP (1) be periodically stationary such that (4) holds.

  * (ii)

Let $\\{\epsilon_{t}\\}_{t\in\mathbb{Z}}$ be a zero-mean periodic white noise
process with period $S$ and positive definite covariance matrices
$\Sigma_{\epsilon}(s),s=1\dots,S$ with

$\epsilon_{Sn+s}=H_{0}(s)w_{Sn+s},\;\;n\in\mathbb{Z},$

where $\\{w_{t}\\}_{t\in\mathbb{Z}}$ is a strictly stationary white noise
process with $E(w_{t}w_{t}^{\prime})=\bm{I}_{m}$, $t\in\mathbb{Z}$, and non-
singular matrix $H_{0}(s)$ satisfying
$\Sigma_{\epsilon}(s)=H_{0}(s)H_{0}(s)^{\prime}$, $s=1,\dots,S$.

  * ($iii)$

The process $\\{w_{t}\\}_{t\in\mathbb{Z}}$ is $\alpha$-mixing with
$\sum_{h=1}^{\infty}h^{\rho-2}(\alpha(h))^{\delta/(2\rho-2+\delta)}<\infty$
for $\rho=4$ and some $\delta>0$, where
$\alpha(h)=\sup_{A\in\mathcal{F}_{-\infty}^{0},B\in\mathcal{F}_{h}^{\infty}}|P(A\cap
B)-P(A)P(B)|$, $h\in\mathbb{N}$ denote the $\alpha$-mixing coefficients of
$\\{w_{t}\\}_{t\in\mathbb{Z}}$ with
$\mathcal{F}_{-\infty}^{0}=\sigma(\dots,w_{-2},w_{-1},w_{0})$ and
$\mathcal{F}_{h}^{\infty}=\sigma(w_{h},w_{h+1},\dots)$. Further, we assume
that $E|w_{t}|^{\rho+2\delta}_{\rho+2\delta}<\infty$, where $|.|_{p}$ denotes
the entry-wise $p$-norm.

Instead of the common i.i.d. assumption for the white noise process
$\\{w_{t}\\}_{t\in\mathbb{Z}}$, we assume the less restrictive $\alpha$-mixing
assumption. Deriving theory under this more general assumption of innovations
that are uncorrelated, but possibly dependent, is important in practice. For
instance, as demonstrated in our real-data application in Figure 5 in Section
6 below, the (P)VAR residuals are likely to be non-linearly dependent.
Following Jentsch & Lunsford (2022), we impose a typical $\alpha$-mixing
assumption, which implies absolute summability of (joint) cumulants up to
order $4$. Periodic weak white noise processes cover a large class of
uncorrelated, but possibly dependent periodically (strictly) stationary
processes and allow, e.g., for periodic conditional heteroskedasticity, see
Francq et al. (2011), Bibi & Ghezal (2016) and Boubacar Maïnassara & Ursu
(2023). While Boubacar Maïnassara & Ursu (2023) already showed asymptotic
normality in the weak PVAR case for linearly restricted least squares
estimators _within_ a given season, we generalize this result for linearly
constrained least squares estimators that also allow for restrictions _across_
seasons. Additionally, we provide asymptotic results for covariance estimators
. Under Assumption 3.1, we obtain a joint CLT of $\widehat{\beta}_{res}$ and
$\widehat{\sigma}$.

######  Theorem 3.1 (Joint CLT with Periodic Weak White Noise).

Under Assumption 3.1, we have

$\sqrt{T}\begin{pmatrix}\widehat{\beta}_{res}-\beta\\\
\widehat{\sigma}-\sigma\end{pmatrix}\overset{d}{\to}\mathcal{N}(\bm{0},SV),\;\;\;V=\begin{pmatrix}V^{(1,1)}&V^{(1,2)}\\\
V^{(2,1)}&V^{(2,2)}\end{pmatrix},$

where the submatrices $V^{(1,1)},V^{(2,1)},V^{(1,2)}$ and $V^{(2,2)}$ are of
dimension $(m\sum_{s=1}^{S}(mp(s)+1)\times
m\sum_{s=1}^{S}(mp(s)+1)),(S\widetilde{m}\times
m\sum_{s=1}^{S}(mp(s)+1)),(m\sum_{s=1}^{S}(mp(s)+1)\times S\widetilde{m})$ and
$(S\widetilde{m}\times S\widetilde{m})$ and are provided in supplementary
material A.

Due to the seasonal structure of the estimators $\widehat{\beta}_{res}$ and
$\widehat{\sigma}$, the submatrices of $V$ are composed of further season-wise
submatrices leading to a tedious and lengthy notation. In particular, the
submatrices of $V$ depend on infinite sums of $\tau_{a,b,c}(s_{1},s_{2})$ and
$\kappa_{a,b}(s_{1},s_{2})$, $s_{1},s_{2}=1,\dots,S$ and $a,b,c\in\mathbb{Z}$
given in (11) and (12), which are cumbersome to estimate, and the submatrix
$V^{(2,1)}=V^{(1,2)\prime}$ is generally not the zero-matrix such that the
estimators $\widehat{\beta}_{res}$ and $\widehat{\sigma}$ become
asymptotically dependent. We recommend bootstrap methods to be discussed in
Section 5 to approximate the limiting variances when conducting inference. If
we impose stronger assumptions such as e.g. martingale difference sequences or
i.i.d. assumptions, the submatrices of $V$ simplify a lot. The i.i.d. case is
stated in Remark 3.3 below.

The CLT in Theorem 3.1 can be either represented with rate of convergence
$\sqrt{T}$ and $S$ inflating the limiting variance or with rate of convergence
$\sqrt{N}$ and without factor $S$ in the limiting variance. Hence, compared to
the corresponding VAR case, we have to pay a price that can be expressed
either in terms of convergence rate or uncertainty. This is not surprising at
all as the ability of PVARs to capture seasonalities is associated with larger
parameterizations than in commonly used VARs. However, for the special case of
$S=1$, where the weak PVAR collapses to a weak VAR, we obtain the same results
as in Francq & Raïssi (2007) and Brüggemann et al. (2016).

######  Remark 3.1.

Under Assumption 3.1 and if no restrictions across seasons are imposed, that
is, we can write the linear restrictions as $\beta(s)=R(s)\gamma(s)+r(s)$,
$s=1,\ldots,S$, where $R(s)$ is of dimension $(m^{2}p(s)+1)\times{M(s)}$ and
$r(s)$ of dimension $m^{2}p(s)+1$ for $s=1,\dots,S$, our asymptotic results
cover the asymptotic results stated in Theorem 3.1 of Boubacar Maïnassara &
Ursu (2023) as a special case.

######  Remark 3.2.

Let Assumption 3.1 hold with $H_{0}(s)=H$ and $p(s)=p$, $s=1,\dots,S$, and let
$R=\bm{1}_{S}\otimes\bm{I}_{m(mp+1)}$ and $r=\bm{0}_{Sm(mp+1)}$. Then, we get
the results on standard VARs for weak white noise assumptions used by
Brüggemann et al. (2016) and Jentsch & Lunsford (2022).

Remark 3.1 states that in the periodic weak white noise case, we get the same
asymptotic results as Boubacar Maïnassara & Ursu (2023), if we do not impose
restrictions across seasons. However, in contrast to Boubacar Maïnassara &
Ursu (2023), we do not obtain block diagonality of $V^{(1,1)}$. Remark 3.2
states that Theorem 3.1 covers also the VAR framework, if all periodic VAR
coefficients are restricted to be equal across all seasons.

######  Remark 3.3.

Now, we consider the special case of a strong PVAR under general linear
constraints, which is also covered by Theorem 3.1. For this purpose, we impose
a periodic strong white noise (Ursu & Duchesne, 2009) assumption on
$\\{y_{t}\\}_{t\in\mathbb{Z}}$ achieved by replacing $(ii)$ and $(iii)$ in
Assumption 3.1 by $(ii)^{\prime}$ and $(iii)^{\prime}$, where

  * $(ii)^{\prime}$

Let $\\{\epsilon_{t}\\}_{t\in\mathbb{Z}}$ be a zero-mean periodic white noise
process with period $S$ and positive definite covariance matrices
$\Sigma_{\epsilon}(s),s=1\dots,S$, and $\\{\epsilon_{t}\\}_{t\in\mathbb{Z}}$
can be represented as

$\epsilon_{Sn+s}=H_{0}(s)w_{Sn+s},\quad n\in\mathbb{Z},$

where $\\{w_{t}\\}_{t\in\mathbb{Z}}$ with
$w_{t}=(w_{t,1},\dots,w_{t,m})^{\prime}$ is an i.i.d. white noise with
$E(w_{t}w_{t}^{\prime})=\bm{I}_{m}$, $t\in\mathbb{Z}$ and the non-singular
matrix $H_{0}(s)$ satisfies $\Sigma_{\epsilon}(s)=H_{0}(s)H_{0}(s)^{\prime}$,
$s=1,\dots,S$.

  * $(iii)^{\prime}$

Let $E|w_{t}|^{4}_{4}<\infty$, where $|.|_{p}$ denotes the entry-wise
$p$-norm.

In this case, we obtain a joint CLT of the PVAR estimators
$\widehat{\beta}_{res}$ and $\widehat{\sigma}$ with simplified limiting
variances. Details can be found in Corollary B.1 in the supplementary
material.

##  4 Structural PVAR Analysis

In this section, we derive the structural-form PVAR model (SPVAR) and derive
asymptotic properties of (structural) impulse responses. As in the VAR setup,
PVAR models do not capture contemporaneous effects between economic variables.
This structure remains in the error term as (possibly) periodic
contemporaneous correlation. Structural models are able to capture these
effects by transforming the reduced-form error terms into orthogonal
structural shocks. Thus, we impose the following assumptions on the variance-
covariance matrices of the periodic white noise error term

$\displaystyle\Sigma_{\epsilon}(s)=H_{0}(s)H_{0}(s)^{\prime},\quad
s=1,\dots,S,$ (13)

where $H_{0}(s)\in\mathbb{R}^{m\times m}$ is invertible due to the non-
singularity of $\Sigma_{\epsilon}(s)$. Then, transforming the reduced-form
PVAR model by left-multiplying with $H_{0}^{-1}(s)$ yields the SPVAR model

$\displaystyle H_{0}^{-1}(s)y_{Sn+s}$
$\displaystyle=H_{0}^{-1}(s)A_{1}(s)y_{Sn+s-1}+\dots+H_{0}^{-1}(s)A_{p(s)}(s)y_{Sn+s-p(s)}+w_{Sn+s},$

where $w_{Sn+s}=H_{0}^{-1}(s)\epsilon_{Sn+s}$ denotes the structural error
term satisfying $Var(w_{Sn+s})=Var(H_{0}^{-1}(s)\epsilon_{Sn+s})=\bm{I}_{m}$,
$s=1,\dots,S$.

As the variance-covariance matrix of the structural shocks does not depend on
the season $s$, the structural error term $\\{w_{t}\\}_{t\in\mathbb{Z}}$
becomes a strictly stationary white noise process with the identity matrix as
variance-covariance matrix. This result is fairly consistent with SVAR theory
since it implies that the endogenous variables are driven by shocks that are
white noise. Moreover, in line with Assumption 3.1, this result implies that
the whole periodic structure of the variables is captured by the SPVAR model.

Structural shocks are often identified based on economic theories and in order
to make these shocks economically interpretable. Common SVAR identification
techniques as short-run, long-run, proxy, sign and inequality restrictions,
see Kilian & Lütkepohl (2017) for an overview, can also be used for PVARs to
identify economically interpretable structural shocks222Compared to the SVAR
setup, (13) is uniquely solvable, if at least $m(m-1)/2$ restrictions are
imposed on
${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}H}_{0}(s)^{-1}$
for each season $s=1,\dots,S$. A trivial way to identify structural shocks is
defining
${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}H}_{0}(s)^{-1}$
as the lower-triangular Cholesky factor of $\Sigma_{\epsilon}(s)$,
$s=1,\dots,S$. However, this recursive short-run identification method does
not necessarily lead to economically interpretable structural shocks. . As
these identification methods have to be imposed on _each_ season $s=1,\dots
S$, this can cause biased results since applying structural identification
techniques on the PVAR can be much more restrictive than the same
identification scheme used on a VAR.

Structural impulse responses at lag $k$ are defined as a reaction of
$y_{Sn+s+k}$ in response to a one-time impulse in $w_{Sn+s}$. Hence, following
Kilian & Lütkepohl (2017), structural PVAR impulse responses for each season
$s=1,\dots,S$ are defined as

$\Theta^{SIR}_{k}(s)=\frac{\partial y_{Sn+s+k}}{\partial w_{Sn+s}},\quad
k\in\mathbb{N}_{0}.$

For $i,j=1,\dots,m$, the element in the $i$-th row and $j$-th column of
$\Theta^{SIR}_{k}(s)$ denotes the structural impulse response of the $i$-th
component of $y_{Sn+s+k}$ to a unit shock in the $j$-th component of
$w_{Sn+s}$ at season $s$. The structural impulse responses can also be deduced
by transforming the moving average representation of $y_{Sn+s}$ as follows

$\displaystyle
y_{Sn+s}=\mu(s)+\sum_{k=0}^{\infty}\Phi_{k}(s)\epsilon_{Sn+s-k}=\mu(s)+\sum_{k=0}^{\infty}\Theta_{k}(s)w_{Sn+s-k},$

where $\Theta_{k}(s)=\Phi_{k}(s)H_{0}(s-k)$ are the structural moving average
coefficient matrices indicating the effect of $w_{Sn+s-k}$ on $y_{Sn+s}$. Note
that $H_{0}(s)=H_{0}(Sn+s)$ for all $n\in\mathbb{Z}$ and $s=1,\dots,S$. The
structural impulse responses $\Theta^{SIR}_{k}(s)$ can expressed analytically
as

$\Theta^{SIR}_{k}(s)=\Theta_{k}(s+k)=\Phi_{k}(s+k)H_{0}(s)=\Phi^{IR}_{k}(s)H_{0}(s).$

In writing down the asymptotic properties of the impulse responses
$\Phi_{k}^{IR}(s)$ and the structural impulse responses $\Theta_{k}^{SIR}(s)$,
$s=1,\dots,S$, we use the higher dimensional VAR representation of the PVAR in
(3). The VAR impulse responses $\Pi_{h}^{IR}$, $h\in\mathbb{N}_{0}$, in (6),
consist of the PVAR impulse responses $\Phi_{k}^{IR}(s)$ for all $k,s$
satisfying $hS<k+s\leq(h+1)S$, while the structural impulse responses of the
higher dimensional VAR denoted by $\Psi_{h}^{SIR}$, $h\in\mathbb{N}_{0}$
consist of $\Theta_{k}^{SIR}(s)$ for all $k,s$ that satisfy
$hS<k+s\leq(h+1)S$. Precisely, we have

$\Psi_{h}^{SIR}=\Pi_{h}^{IR}H_{0}=\begin{pmatrix}\Theta_{Sh}^{SIR}(1)&\Theta_{Sh-1}^{SIR}(2)&\dots&\Theta_{Sh-S+1}^{SIR}(S)\\\
\Theta_{Sh+1}^{SIR}(1)&\Theta_{Sh}^{SIR}(2)&\ddots&\vdots\\\
\vdots&\ddots&\ddots&\Theta_{Sh-1}^{SIR}(S)\\\
\Theta_{Sh+S-1}^{SIR}(1)&\dots&\Theta_{Sh+1}^{SIR}(S-1)&\Theta_{Sh}^{SIR}(S)\end{pmatrix},\;\;\;h\in\mathbb{N}_{0},$

where
$H_{0}=\text{diag}\bigl{(}H_{0}(1),H_{0}(2),\dots,H_{0}(S)\bigr{)}\in\mathbb{R}^{Sm\times
Sm}$. Let $\widehat{\Pi}_{h}^{IR}$ and $\widehat{\Psi}_{h}^{SIR}$ be the
empirical counterparts of $\Pi_{h}^{IR}$ and $\Psi_{h}^{SIR}$ for
$h\in\mathbb{N}_{0}$, respectively. Note that $\Pi_{h}^{IR}$,
$h\in\mathbb{N}_{0}$ can be represented as a continuously differentiable
function of the PVAR coefficients $\beta$, while the structural impulse
responses $\Psi_{h}^{SIR}$ can be represented as a continuously differentiable
function depending on $\beta$ and $\sigma$. Consequently, from Theorem 3.1 and
following Lütkepohl (2005), Proposition 3.6, we get immediately the following
result.

######  Theorem 4.1 (CLT of Impulse Responses).

Under Assumption 3.1, we have

$\displaystyle\sqrt{T}\bigl{(}vec\\{\widehat{\Pi}_{h}^{IR}-\Pi_{h}^{IR}\\}\bigr{)}\overset{d}{\to}\mathcal{N}\bigl{(}\bm{0},S\;\Sigma_{h}^{\Pi}\bigr{)},$
$\displaystyle\sqrt{T}\bigl{(}vec\\{\widehat{\Psi}_{h}^{SIR}-\Psi_{h}^{SIR}\\}\bigr{)}\overset{d}{\to}\mathcal{N}\bigl{(}\bm{0},S\;\Sigma_{h}^{\Psi}\bigr{)},\;\;\;h\in\mathbb{N}_{0},$

where

$\displaystyle\Sigma_{h}^{\Pi}=G^{\Pi}_{h}V^{(1,1)}G_{h}^{\Pi\prime}\quad\text{and}\quad\Sigma_{h}^{\Psi}=F_{h}V^{(1,1)}F_{h}^{\prime}+D_{h}V^{(2,2)}D_{h}^{\prime}+F_{h}V^{(2,1)\prime}D_{h}^{\prime}+D_{h}V^{(2,1)}F_{h}^{\prime}$

with $G_{h}=\frac{\partial vec\\{\Pi_{h}^{IR}\\}}{\partial\beta^{\prime}}$,
$F_{h}=\frac{\partial vec\\{\Psi_{h}^{SIR}\\}}{\partial\beta^{\prime}}$ and
$D_{h}=\frac{\partial vec\\{\Psi_{h}^{SIR}\\}}{\partial\sigma^{\prime}}$.

The closed-form solutions of the derivatives $G_{h},F_{h}$ and $D_{h}$ are
stated in supplemental appendix C. In the strong PVAR case, the same
statements as in Theorem 4.1 hold with the matrices $V^{(1,1)},V^{(2,1)}$ and
$V^{(2,2)}$ replaced by $V^{(1,1)}_{iid},V^{(2,1)}_{iid}$ and
$V^{(2,2)}_{iid}$, respectively.

##  5 Bootstrap Inference in PVARs

Analytical expressions of the limiting variances of the PVAR estimators and
the (structural) impulse responses are very complex and cumbersome to
estimate. Therefore, we propose residual-based bootstrap methods that are
suitable for periodically correlated processes to approximate the (limiting)
distributions of the estimators and to construct confidence intervals for
(structural) impulse response. Further, we prove their bootstrap consistency.

###  5.1 Bootstrap Schemes

In this section, we discuss different residual-based bootstrap methods that
are suitable for conducting inference in different PVAR setups. In Section
5.1.1, we describe residual-based bootstrap schemes that are tailored for the
assumption of a periodic white noise, while Section 5.1.2 contains simplified
versions that are sufficient for periodic strong white noise.

####  5.1.1 Residual-based seasonal block bootstrapping

We describe residual-based seasonal block bootstrap methods for PVAR processes
that work under Assumption 3.1, where the PVAR error term process
$\\{\epsilon_{t}\\}_{t\in\mathbb{Z}}$ is assumed to be periodic, possibly
_weak_ white noise. In this context, block bootstrapping is required to mimic
the potential non-linear serial dependence, while a seasonal version is
required because $\\{\epsilon_{t}\\}_{t\in\mathbb{Z}}$ is allowed to be
periodically correlated. Based on sample values $y_{1},\dots,y_{T}$ with
$T=SN$ and pre-sample values $y_{s^{\prime}},\dots,y_{0}$, suppose that
estimates of the PVAR($\bm{p})$ parameters $\beta$ and $\sigma$ are given
under general linear constraints $\beta=R\gamma+r$ as stated in Section 3.1.

The algorithm is as follows. We initialize the algorithm by choosing a
positive integer block size $b<T$ and define $l=\lceil T/b\rceil$, where
$\lceil.\rceil$ denotes the ceiling function. Next, we define $(m\times b)$
blocks
$\widehat{\mathcal{E}}_{i}=(\widehat{\epsilon}_{i},\dots,\widehat{\epsilon}_{i+b-1})$
for $i=1,\dots,T-b+1$. Then, we proceed as follows:

  * $(1)$

For $t=1,b+1,2b+1,\dots,(l-1)b+1$, define the block of pseudo residuals by
$\widehat{\mathcal{E}}^{*}_{t}=\widehat{\mathcal{E}}_{k_{t}}$, where
$k_{1},k_{b+1},k_{2b+1},\dots$ are independent, discrete uniform random
variables on the set

$\\{t-SR_{1,T},t-S(R_{1,T}-1),\dots,t-S,t,t+S,\dots,t+S(R_{2,n}-1),t+SR_{2,n}\\}$

with $R_{1,T}=\lfloor\frac{t-1}{S}\rfloor$ and
$R_{2,T}=\lfloor\frac{T-b-t}{S}\rfloor$, where $\lfloor.\rfloor$ denotes the
floor function.

  * $(2)$

Collect the blocks
$(\widehat{\mathcal{E}}^{*}_{1},\dots,\widehat{\mathcal{E}}^{*}_{(l-1)b+1})$
and remove the last $lb-T$ observations, leaving a sample
$\widehat{\epsilon}_{1}^{*},\dots,\widehat{\epsilon}_{T}^{*}$ of size $T$.

  * $(3)$

Set $(y^{*}_{s^{\prime}},\dots,y^{*}_{0})=(y_{s^{\prime}},\dots,y_{0})$ and
use estimates of $\beta$ and $\sigma$ and
$\widehat{\epsilon}_{1}^{*},\dots,\widehat{\epsilon}_{T}^{*}$ to recursively
compute bootstrap observations $y^{*}_{1},\dots,y^{*}_{T}$ according to

$\displaystyle
y^{*}_{Sn+s}=\widehat{\nu}(s)+\widehat{A}_{1}(s)y^{*}_{Sn+s-1}+\dots+\widehat{A}_{p(s)}(s)y^{*}_{Sn+s-p(s)}+\widehat{\epsilon}^{*}_{Sn+s}$

for $s=1,\dots,S$ and $n=0,\dots,N-1$.

  * $(4)$

Reestimate $\beta$, $\sigma$, $\Pi_{h}^{IR}$ and $\Psi_{h}^{SIR}$ to get
$\widehat{\beta}_{res}^{*}$, $\widehat{\sigma}^{*}$, $\widehat{\Pi}_{h}^{IR*}$
and $\widehat{\Psi}_{h}^{SIR*}$ from the bootstrap (pre)sample
$(y^{*}_{s^{\prime}},\ldots,y_{0}^{*}),y_{1}^{*}\dots,y^{*}_{T}$ and use the
same linear constraints and identification scheme as in the calculation of the
estimates $\widehat{\beta}_{res}$, $\widehat{\sigma}$,
$\widehat{\Pi}_{h}^{IR}$ and $\widehat{\Psi}_{h}^{SIR}$.

In the bootstrap algorithm, a seasonal block bootstrap proposed by Dudek et
al. (2014, 2016) is applied directly to the PVAR residuals. Alternatively, a
non-seasonal resampling variant can also be used. To initialize this
bootstrap, we choose a positive integer block size $b<T$ and define $l=\lceil
T/b\rceil$. Next, we define $(m\times b)$ blocks
$\widehat{\mathcal{U}}_{i}=(\widehat{u}_{i},\dots,\widehat{u}_{i+b-1})$,
$i=1,\dots,T-b+1$, where
$\widehat{u}_{Sn+s}=\widehat{\Sigma}_{\epsilon}^{-1/2}(s)\widehat{\epsilon}_{Sn+s}$
denote the _standardized_ residuals for $s=1,\dots,S$ and $n=0,\dots,N-1$.
Then, the non-seasonal variant is obtained by replacing steps $(1)-(2)$ by
$(1^{\prime})-(2^{\prime})$, where

  * $(1^{\prime})$

For $t=1,\dots,l$, define the block of pseudo standardized residuals by
$\widehat{\mathcal{U}}^{*}_{t}=\widehat{\mathcal{U}}_{k_{t}}$, where
$k_{1},k_{2},\dots$ are independent, discrete uniform random variables on
$\\{1,\dots,T-b+1\\}$.

  * $(2^{\prime})$

Collect the blocks
$(\widehat{\mathcal{U}}^{*}_{1},\dots,\widehat{\mathcal{U}}^{*}_{l})$ and
remove the last $lb-T$ observations, leaving a sample
$\widehat{u}_{1}^{*},\dots,\widehat{u}_{T}^{*}$ of size $T$ and generate
bootstrap PVAR residuals
$\widehat{\epsilon}_{1}^{*},\dots,\widehat{\epsilon}_{T}^{*}$ by

$\widehat{\epsilon}_{Sn+s}^{*}=\widehat{\Sigma}_{\epsilon}^{1/2}(s)\widehat{u}_{Sn+s}^{*},\;\;s=1,\dots,S,\;n=0,\dots,N-1.$

In this alternative bootstrap method, the periodic residuals are first
standardized and then a classical moving-block bootstrap approach, as
described, e.g., in Jentsch & Lunsford (2022), is applied to the adequately
standardized residuals. As both bootstrap algorithms need a suitably chosen
block length $b$ as initialization, we refer to Nordman (2009) and Bertail &
Dudek (2022), who developed optimal choices for the block length of non-
seasonal and seasonal block bootstrap methods, respectively. With $b=1$, both
block bootstrap methods collapse to the residual-based independent bootstraps
discussed in Section 5.1.2.

By repeating the bootstrap schemes $L$ times, where $L$ is large, we get
estimates of the standard errors of the PVAR parameters and confidence
intervals of the structural impulse responses. We construct the confidence
intervals of size $\alpha$ of the structural impulse responses
$\Psi_{h}^{SIR}$, $h\in\mathbb{Z}$ are produced according to

$\displaystyle
CI_{\alpha}=[\widehat{\Psi}_{h}^{SIR}+s^{*}_{\alpha/2,h}-s^{*}_{0.5,h},\widehat{\Psi}_{h}^{SIR}+s^{*}_{1-\alpha/2,h}-s^{*}_{0.5,h}],$
(14)

where $s^{*}_{\alpha,h}$ denotes the empirical $\alpha$-quantile of the $L$
bootstrap values
$\widehat{\Psi}_{h}^{SIR*(1)},\ldots,\widehat{\Psi}_{h}^{SIR*(L)}$.
Alternatively, standard percentile, Hall’s percentile or Hall’s studentized
intervals can also be used. However, our simulation study shows that the
method for the construction of confidence intervals proposed in (14) provides
slightly more accurate coverage.

The intuition why for periodic weak white noise errors a block bootstrap
approach is needed is that an independent bootstrap is not able to capture
possibly present (non-linear) serial dependencies of the weak white noise
process. Note that following Theorem 3.1, the limiting variances $V$ consist
of infinite sums of $\tau_{a,b,c}(s_{1},s_{2})$ and
$\kappa_{a,b}(s_{1},s_{2})$, $s_{1},s_{2}=1,\dots,S$ and $a,b,c\in\mathbb{Z}$
and the independent drawing of the residuals destroys this serial dependence.
Since residual-based block bootstrap approaches resample joint blocks of the
residuals, they capture asymptotically the entire serial dependence structure
and are therefore able to mimic the limiting variances $V$ correctly.

####  5.1.2 Residual-based periodic independent bootstrapping

If Assumption 3.1 with $(ii)$ and $(iii)$ replaced by $(ii)^{\prime}$ and
$(iii)^{\prime}$ holds, where the PVAR innovation process
$\\{\epsilon_{t}\\}_{t\in\mathbb{Z}}$ is assumed to be periodic strong white
noise, simplified residual-based independent bootstrap techniques can be used.
In Section 5.1.1, by setting $b=1$ in steps $(1)-(2)$, the seasonal version of
the residual-based independent bootstrap is obtained, while the non-seasonal
version is obtained by setting $b=1$ in steps $(1^{\prime})-(2^{\prime})$.

The residual-based independent bootstrap methods provide valid confidence
intervals only when Assumption 3.1 with $(ii)$ and $(iii)$ replaced by
$(ii)^{\prime}$ and $(iii)^{\prime}$ hold. However, as illustrated in Figure
5, the periodic strong white noise assumptions may be violated in
applications. Hence, we recommend to use residual-based block bootstrap
methods for inference of periodically correlated time series.

###  5.2 Bootstrap Theory for PVARs

We show bootstrap consistency for the residual-based seasonal block bootstrap
methods discussed in Sections 5.1.1 and 5.1.2 under periodic weak and strong
white noise assumptions, respectively. Addressing first the weak PVAR case, we
need the following assumption.

######  Assumption 5.1 (Periodic Weak White Noise with finite eighth-order
Moments).

The weak white noise process $\\{w_{t}\\}_{t\in\mathbb{Z}}$ as defined in
Assumption 3.1 (ii) is $\alpha$-mixing such that Assumption 3.1 (iii) holds
for $\rho=8$.

For weak white noise processes, such an assumption is made, e.g., in Gonçalves
& Kilian (2007) to guranatee summability of cumulants up to order eight to be
able to prove consistency of wild and pairwise bootstrap methods for
AR($\infty$) processes, in Brüggemann et al. (2016) to prove consistency of a
residual-based moving block bootstrap for VAR($p$) models, and in Jentsch &
Lunsford (2019, 2022) to prove consistency of a residual-based moving-block
bootstrap for Proxy SVARs. We now show that the residual-based seasonal block
bootstrap can approximate the limiting distributions of
$\sqrt{T}\bigl{(}\bigl{(}\widehat{\beta}_{res}-\beta\bigr{)}^{\prime},\bigl{(}\widehat{\sigma}-\sigma\bigr{)}^{\prime}\bigr{)}^{\prime}$,
$\sqrt{T}\bigl{(}vec\\{\widehat{\Pi}_{h}^{IR}-\Pi_{h}^{IR}\\}\bigr{)},$ and
$\sqrt{T}\bigl{(}vec\\{\widehat{\Psi}_{h}^{SIR}-\Psi_{h}^{SIR}\\}\bigr{)},$
$h\in\mathbb{N}_{0}$ derived in Theorems 3.1 and 4.1.

######  Theorem 5.1 (Residual-Based Seasonal Block Bootstrap Consistency).

Suppose Assumptions 3.1 and 5.1 hold and the residual-based seasonal block
bootstrap schemes discussed in Section 5.1.1 are used to obtain bootstrap
estimators
$\widehat{\beta}^{*}_{res},\widehat{\sigma}^{*},\widehat{\Pi}_{h}^{IR*}$ and
$\widehat{\Psi}_{h}^{SIR*}$ for $h\in\mathbb{N}_{0}.$ If $b\to\infty$ such
that $b^{3}/T\to 0$ as $T\to\infty$, we get

$\displaystyle\sup_{x\in\mathbb{R}^{m_{1}}}\Big{|}P^{*}\Bigl{(}\sqrt{T}\bigl{(}\bigl{(}\widehat{\beta}^{*}_{res}-\widehat{\beta}_{res}\bigr{)}^{\prime},\bigl{(}\widehat{\sigma}^{*}-\widehat{\sigma}\bigr{)}^{\prime}\bigr{)}^{\prime}\leq
x\Bigr{)}-P\Bigl{(}\sqrt{T}\bigl{(}\bigl{(}\widehat{\beta}_{res}-\beta\bigr{)}^{\prime},\bigl{(}\widehat{\sigma}-\sigma\bigr{)}^{\prime}\Bigr{)}^{\prime}\leq
x\bigr{)}\Big{|}\overset{p}{\to}0,$

where $m_{1}=m\sum_{s=1}^{S}(mp(s)+1)+S\widetilde{m}$, and, for
$h\in\mathbb{Z}$, we have

$\displaystyle\sup_{x\in\mathbb{R}^{(Sm)^{2}}}\Big{|}P^{*}\Bigl{(}\sqrt{T}\bigl{(}vec\\{\widehat{\Pi}_{h}^{IR*}-\widehat{\Pi}_{h}^{IR}\\}\bigr{)}\leq
x\Bigr{)}-P\Bigl{(}\sqrt{T}\bigl{(}vec\\{\widehat{\Pi}_{h}^{IR}-\Pi_{h}^{IR}\\}\bigr{)}\leq
x\Bigr{)}\Big{|}\overset{p}{\to}0,$
$\displaystyle\sup_{x\in\mathbb{R}^{(Sm)^{2}}}\Big{|}P^{*}\Bigl{(}\sqrt{T}\bigl{(}vec\\{\widehat{\Psi}_{h}^{SIR*}-\widehat{\Psi}_{h}^{SIR}\\}\bigr{)}\leq
x\Bigr{)}-P\Bigl{(}\sqrt{T}\bigl{(}vec\\{\widehat{\Psi}_{h}^{SIR}-\Psi_{h}^{SIR}\\}\bigr{)}\leq
x\Bigr{)}\Big{|}\overset{p}{\to}0,$

where $P^{*}$ denotes the probability measure induced by the residual-based
seasonal block bootstrap methods. Note that $x\leq y$ for some
$x,y\in\mathbb{R}^{K}$ stands for $x_{i}\leq y_{i}$ for all $i=1,\dots,K$.

The proof of Theorem 5.1 can be found in supplemental appendix D. For the
consistency result of the (structural) impulse responses, we use that the
impulse responses $\Pi_{h}^{IR}$ and structural impulse responses
$\Psi_{h}^{SIR}$, $h\in\mathbb{Z},$ can be represented as continuously
differentiable functions of $\beta$ and $\sigma$. It is worth noting that in
the case, where the more restrictive periodic strong white noise assumption
applies, the independent bootstrap, i.e., the block bootstrap with block
length $b=1$, is sufficient.

######  Corollary 5.1 (Residual-Based independent Bootstrap Consistency).

Suppose Assumption 3.1 with $(ii)$ and $(iii)$ replaced by $(ii)^{\prime}$ and
$(iii)^{\prime}$ and Assumption 5.1 hold333Note that in the periodic strong
white noise case, Assumption 5.1 boils down to the condition
$E|w_{t}|^{8}_{8}<\infty$. and suppose that the residual-based independent
bootstrap schemes discussed in Section 5.1.2 are used to obtain bootstrap
estimators
$\widehat{\beta}^{*}_{res},\widehat{\sigma}^{*},\widehat{\Pi}_{h}^{IR*}$ and
$\widehat{\Psi}_{h}^{SIR*}$ for $h\in\mathbb{N}_{0}$. Then, as $T\to\infty$,
we get

$\displaystyle\sup_{x\in\mathbb{R}^{m_{1}}}\Big{|}P^{*}\Bigl{(}\sqrt{T}\bigl{(}\bigl{(}\widehat{\beta}^{*}_{res}-\widehat{\beta}_{res}\bigr{)}^{\prime},\bigl{(}\widehat{\sigma}^{*}-\widehat{\sigma}\bigr{)}^{\prime}\bigr{)}^{\prime}\leq
x\Bigr{)}-P\Bigl{(}\sqrt{T}\bigl{(}\bigl{(}\widehat{\beta}_{res}-\beta\bigr{)}^{\prime},\bigl{(}\widehat{\sigma}-\sigma\bigr{)}^{\prime}\Bigr{)}^{\prime}\leq
x\bigr{)}\Big{|}\overset{p}{\to}0,$

where $m_{1}=m\sum_{s=1}^{S}(mp(s)+1)+S\widetilde{m}$, and, for
$h\in\mathbb{Z}$, we have

$\displaystyle\sup_{x\in\mathbb{R}^{(Sm)^{2}}}\Big{|}P^{*}\Bigl{(}\sqrt{T}\bigl{(}vec\\{\widehat{\Pi}_{h}^{IR*}-\widehat{\Pi}_{h}^{IR}\\}\bigr{)}\leq
x\Bigr{)}-P\Bigl{(}\sqrt{T}\bigl{(}vec\\{\widehat{\Pi}_{h}^{IR}-\Pi_{h}^{IR}\\}\bigr{)}\leq
x\bigr{)}\Big{|}\overset{p}{\to}0,$
$\displaystyle\sup_{x\in\mathbb{R}^{(Sm)^{2}}}\Big{|}P^{*}\Bigl{(}\sqrt{T}\bigl{(}vec\\{\widehat{\Psi}_{h}^{SIR*}-\widehat{\Psi}_{h}^{SIR}\\}\bigr{)}\leq
x\Bigr{)}-P\Bigl{(}\sqrt{T}\bigl{(}vec\\{\widehat{\Psi}_{h}^{SIR}-\Psi_{h}^{SIR}\\}\bigr{)}\leq
x\bigr{)}\Big{|}\overset{p}{\to}0,$

where $P^{*}$ denotes the probability measure induced by the residual-based
independent bootstrap.

##  6 Real Data Application

In this section, we compare the results of a structural PVAR applied to not
seasonally adjusted macroeconomic data with the corresponding results of a
structural VAR applied to corresponding seasonally adjusted data. The
variables included in the model are monthly industrial production (IP) as a
measure of real economic activity, CPI inflation (INF) and federal funds rate
(FFR) from 1968-2019 taken from the FRED database. In the macroeconomic
literature, it is common to analyze the above-mentioned variables or strongly
related variables, e.g., Bernanke et al. (1997), Peersman (2005) and Kilian &
Lütkepohl (2017). In order to avoid unit roots, we take log differences of IP,
whereas INF and FFR remains untouched. For structural PVAR analysis, IP and
INF are considered on a not seasonally adjusted basis, while their seasonally
adjusted counterparts are considered for SVAR analysis. As FFR is non-
periodic, we use it for both analyses.

Due to the monthly frequency, an obvious choice for the number of seasons
would be $S=12$, resulting in $N=52$ cycles (years). The choice $S=12$ is
strongly supported by the estimates of the SDs and ACFs of IP and INF given in
Figure 1. As demonstrated in Figure 1 and 2 in the introductory section, we
find that IP and INF exhibit strong periodic patterns, while FFR has no
periodic structure at all. The periodic structure of IP is dominated mainly by
sub-year periods and cannot be eliminated by seasonal demeaning alone, while
the periodic structure of INF can be explained mainly by seasonal variation in
the mean. Consequently, seasonal demeaning appears to be not sufficient to
remove the entire seasonal sub-year pattern of IP, while in the case of INF
the periodic structure is almost eliminated completely.

###  6.1 Reduced Form PVAR Model

On the seasonally adjusted data basis, we use a VAR model of order $p=9$ as
the benchmark reduced form model. It captures properly the autoregressive
structure of the monthly macroeconomic data, which is in line with the reduced
form VAR model proposed by Peersman (2005), who use a VAR of order 3 for
quarterly data on oil prices, real GDP, CPI inflation and federal funds rate.
Compared to Peersman (2005), our analysis excludes oil prices, since oil
prices are purely non-periodic and not helpful for the identification of the
structural shocks. For structural PVAR analysis, the reduced form model order
of the PVAR is set to $p(s)=9$, $s=1,\dots,12$ as well. Based on Figures 1 and
2, we treat the FFR to be non-periodic within our periodic setup meaning that
the intercept and all effects of all lagged variables on the FFR are
restricted to be non-seasonal. Further, we found that seasonal demeaning
removes most of the periodic structure of INF. Consequently, we allow for
periodically varying intercepts, but restrict the autoregressive coefficients
of INF to be constant across the seasons for all lags. It can be clearly seen
from Figures 1 and 2 that IP is the only variable in our system that contains
periodic structure which goes beyond seasonality in means. Since the leftover
periodicity of IP after seasonal demeaning is dominated by four-month and
semi-annual periods, we allow for periodically varying intercepts and
autoregressive coefficients up to lag 4 and restrict the remaining
coefficients for lags 5 to 9 to be constant across the seasons. In total,
including intercepts, we estimate 234 parameters, a reasonable number for
monthly data collected from 1968 to 2019. For example, Kilian (2009) and
Kilian & Murphy (2012) estimated similarly parameterized models in a
comparable data setting. The reduced form PVAR(9) with period $S=12$, where
$y^{ip},\pi$ and $i$ denote IP, INF and FFR, is as follows

$\displaystyle
y_{12n+s}=\nu(s)+\sum\limits_{i=1}^{4}A_{i}(s)y_{12n+s-i}+\sum\limits_{i=5}^{9}A_{i}y_{12n+s-i}+\epsilon_{12n+s}$
(15)

for $s=1,\dots,12$ and $n=0,\dots,51$, where

$\displaystyle y_{12n+s}=\begin{pmatrix}y^{ip}_{12n+s}\\\ \pi_{12n+s}\\\
i_{12n+s}\end{pmatrix},\;\;\;\nu(s)=\begin{pmatrix}\nu^{y^{ip}}(s)\\\
\nu^{\pi}(s)\\\ \nu^{i}\end{pmatrix}\quad\text{and}\quad
A_{i}(s)=\begin{pmatrix}a_{i,11}(s)&a_{i,12}&a_{i,13}\\\
a_{i,21}(s)&a_{i,22}&a_{i,23}\\\ a_{i,31}(s)&a_{i,32}&a_{i,33}\end{pmatrix}$

for $i=1,2,3,4$. That is, the first column of the autoregressive coefficient
matrices up to lag 4 are allowed to vary over the seasons, while the remaining
autoregressive coefficient matrices $A_{i}$, $i=5,\ldots,9$ are assumed to be
constant.

Figure 3: ACFs of PVAR (as defined in (15)) residuals
$\widehat{\epsilon}_{1},\widehat{\epsilon}_{2},\dots,\widehat{\epsilon}_{T}$
(top panels) and ACFs of residuals obtained by a VAR(9) fit with seasonal
intercept (bottom panels). For both analyses, not seasonally adjusted monthly
data from 1968-2019 are used.

To check the validity of the imposed restrictions, we compare the residuals
resulting from fitting the restricted PVAR(9) from (15) with those of a VAR(9)
with seasonal intercepts. From Figures 3 and 4, we clearly see that the
restrictions imposed in (15) seem to be very plausible. First, the restricted
PVAR model is able to capture the leftover seasonality in IP and INF much
better than the VAR model. Second, these restrictions do not produce spurious
periodicities. Note that in the PVAR setup, the residuals are allowed to be
periodic. Based on Figure 3, one could argue to choose the model order 12
instead of 9. Despite a significant peak in the ACF at lag 12, the model order
should not be selected to be 12 (or higher) as the autoregressive part of the
model then captures part of the periodic structure, which leads to a
distortion of the results.

###  6.2 Identification

In order to identify aggregate supply $(as)$, aggregate demand $(ad)$ and a
monetary policy $(mp)$ shock, we impose a mixture of short- and long-run
restrictions following Blanchard & Quah (1989), Gali (1992) and Gerlach &
Smets (1995). The assumptions rely on a vertical long-run Philips curve
implying that demand and monetary policy shocks do not have an effect on IP in
the long-run, while aggregate supply shocks are allowed to permanently shift
IP. One could argue that demand shocks may also have long-run effects on IP,
but these long-run effects are usually quite small compared to those of supply
shocks and thus are not likely to affect the estimates significantly. Demand
and monetary policy shocks with permanent effects are therefore labeled as
supply shocks. In order to discriminate between aggregate demand and monetary
policy shocks, following Bernanke & Blinder (1992), we impose that monetary
shocks do not have contemporaneous effects on IP.

Figure 4: Estimated SDs of PVAR (as defined in (15)) residuals
$\widehat{\epsilon}_{1},\widehat{\epsilon}_{2},\dots,\widehat{\epsilon}_{T}$
(black) and of residuals obtained by a VAR(9) fit with seasonal
intercept(green). For both analyses, not seasonally adjusted monthly data from
1968-2019 are used.

###  6.3 Results

In the following, seasonal structural impulse responses using the restricted
SPVAR(9) on not seasonally adjusted data and structural impulse responses
using a non-restricted SVAR(9) on seasonally adjusted data are plotted, in
which the structural shocks are identified using a combination of short- and
long-run restrictions, respectively. The black bold line indicates the
structural impulse response of the SVAR on seasonally adjusted data which can
also be used as the benchmark, while each of the red lines indicates a
structural impulse response generated by the SPVAR for a specific season
$s=1,\dots,12$. The black dashed lines give the 68% confidence interval of the
structural impulse response of the SVAR and the red dashed lines the 68%
confidence intervals of the seasonal impulse responses.

For the construction of confidence intervals of the seasonal impulse
responses, the residual-based seasonal block bootstrap described in Section
5.1.1 has to be used, because Figure 5 shows significant autoregressive
structure in the squared structural shocks such that the strong periodic white
noise assumption appear to not hold. The block length of the residual-based
seasonal block bootstrap is set to $b=7$. This is in line with Bertail & Dudek
(2022). The confidence intervals are constructed by the method proposed in
(14). For the construction of the VAR impulse responses, we use a non-seasonal
block bootstrap in order to construct the confidence intervals.

Figure 5: ACFs of structural PVAR shocks
$\widehat{w}_{1},\widehat{w}_{2},\dots,\widehat{w}_{T}$ (top panels) and of
squared structural PVAR shocks
$\widehat{w}^{2}_{1},\widehat{w}^{2}_{2},\dots,\widehat{w}^{2}_{T}$ (bottom
panels).

Figures 6 and 7 present the seasonal structural impulse responses generated by
a monetary policy shock, which is normalized as an unexpected increase in FFR
by 25 basis points, on IP and INF, respectively. Impulse responses of FFR are
given in supplemental appendix E. We clearly do not see a periodic pattern in
the response of the macro variables after an unexpected monetary policy shock
indicating that the timing of unexpected interest rate hikes or cuts by the
FED has no impact on the variables. Further, we find that the SVAR impulse
responses after a monetary policy shock are approximately given by the mean of
the seasonal structural impulse responses. Interestingly, we do not see
pronounced efficiency losses in the seasonal structural impulse responses in
terms of considerably increased confidence intervals compared to the SVAR
impulse responses. Instead, we even notice some efficiency gains in certain
seasons resulting in more accurate statements; see e.g. the December panel in
Figure 7.

Figure 6: Seasonal structural impulse responses of IP after a monetary policy
shock. SPVAR impulse responses are in red, while SVAR impulse responses are in
black and are constant across the seasons. The corresponding month indicates
the time of occurrence of the shock. Figure 7: Seasonal structural impulse
responses of INF after a monetary policy shock. SPVAR impulse responses are in
red, while SVAR impulse responses are in black and are constant across the
seasons. The corresponding month indicates the time of occurrence of the
shock.  Figure 8: Seasonal structural impulse responses of IP after a negative
supply shock. SPVAR impulse responses are in red, while SVAR impulse responses
are in black and are constant across the seasons. The corresponding month
indicates the time of occurrence of the shock. Figure 9: Seasonal structural
impulse responses of INF after a negative supply shock. SPVAR impulse
responses are in red, while SVAR impulse responses are in black and are
constant across the seasons. The corresponding month indicates the time of
occurrence of the shock.

In terms of the periodic pattern, the structural impulse responses of IP
triggered by a negative supply shock illustrated in Figure 8 show
significantly different results in April (upward shift), July (different
shape) and September (downward shift). In the other seasons, there are hardly
any differences compared to the SVAR impulse responses. As for INF in Figure
9, we clearly see a periodic pattern. The impact on INF of a supply shock
occurring in September appears to be considerably different. It is interesting
to see that in many seasonal impulse responses there is a significant long-
term effect, while this is not the case for the SVAR impulse responses. The
impact of a supply shock on FFR and the structural impulse responses caused by
positive demand shocks are illustrated in supplemental appendix E. The results
clearly show that the effects on the variables vary by season, especially for
shocks that occur in July and September. It is also striking that the impact
of shocks occurring in December is generally smaller and has narrower
confidence intervals than the impact of shocks occurring in the other seasons.
The reasons for these periodic patterns cannot be clearly determined, but we
find these patterns again in various PVAR model variants. In particular, in
July and December, the differences may be explained by certain holiday
effects.

###  6.4 Simulation Study

In Monte Carlo simulation, we examine the seasonal impulse response intervals
constructed by the residual-based seasonal bootstrap methods discussed in
Section 5.1. In doing so, we generate seasonal data by the restricted
structural PVAR used to fit the dynamics between IP, INF and FFR. Hence, in
the simulation study, we generate 3-dimensional $y_{1},\dots,y_{T}$ by

$\displaystyle
y_{12n+s}=\nu(s)+\sum\limits_{i=1}^{4}A_{i}(s)y_{12n+s-i}+\sum\limits_{i=5}^{9}A_{i}y_{12n+s-i}+H_{0}^{-1}(s)w_{12n+s}$

for $s=1,\dots,12$ and $n=0,\dots,51$, where we use the estimated restricted
structural PVAR(9) described in (15) in Section 6.1 as the true underlying
DGP.

We generate the structural shocks $\\{w_{t}\\}_{t\in\mathbb{Z}}$ using
different DGPs, where one variant relies on i.i.d. errors and three variants
make use of differently parametrized GARCH error processes to cover both, the
periodic strong and weak white noise cases. Using
$v_{t}=(v_{1,t},v_{2,t},v_{3,t})^{\prime}\sim
i.i.d.\;\mathcal{N}(\bm{0},\bm{I}_{3}),$ we define
$w_{i,t}=\sigma_{i,t}v_{i,t}$ with
$\sigma^{2}_{i,t}=a_{0}+a_{1}w_{i,t}^{2}+b_{1}\sigma^{2}_{i,t-1},i=1,2,3$ and
$a_{0}=1-a_{1}-b_{1}$ implying that the components of the structural shocks
follow univariate independent GARCH(1,1) processes with
$E(w_{i,t})=1,i=1,2,3$. We use the following GARCH(1,1) specifications of
$\\{w_{t}\\}_{t\in\mathbb{Z}}$ for simulations:

  * (G0)

$a_{1}=0$, $b_{1}=0$ (G1) $a_{1}=0.05$, $b_{1}=0.9$ (G2) $a_{1}=0.3$,
$b_{1}=0.6$ (G3) $a_{1}=0.5$, $b_{1}=0$,

where the case G0 represents the periodic strong white noise case, and
(G1)-(G3) represent variants of the periodic weak white noise case. All these
cases are also considered in Brüggemann et al. (2016).

$b$ | $N$ | $s$ | $as\to y^{ip}$ | $as\to\pi$ | $as\to i$ | $mp\to y^{ip}$ | $mp\to\pi$ | $mp\to i$ | $ad\to y^{ip}$ | $ad\to\pi$ | $ad\to i$
---|---|---|---|---|---|---|---|---|---|---|---
|  | 1 | 0.586 | 0.576 | 0.652 | 0.662 | 0.676 | 0.574 | 0.642 | 0.556 | 0.620
|  | 2 | 0.630 | 0.674 | 0.726 | 0.648 | 0.648 | 0.548 | 0.492 | 0.500 | 0.456
|  | 3 | 0.676 | 0.670 | 0.684 | 0.628 | 0.638 | 0.536 | 0.570 | 0.550 | 0.492
|  | 4 | 0.786 | 0.590 | 0.640 | 0.688 | 0.686 | 0.612 | 0.484 | 0.612 | 0.502
|  | 5 | 0.742 | 0.604 | 0.594 | 0.638 | 0.596 | 0.568 | 0.630 | 0.622 | 0.582
1 | 20 | 6 | 0.610 | 0.584 | 0.638 | 0.712 | 0.694 | 0.640 | 0.656 | 0.662 | 0.648
|  | 7 | 0.738 | 0.516 | 0.578 | 0.646 | 0.622 | 0.608 | 0.572 | 0.630 | 0.562
|  | 8 | 0.584 | 0.638 | 0.658 | 0.642 | 0.666 | 0.604 | 0.512 | 0.458 | 0.468
|  | 9 | 0.542 | 0.598 | 0.600 | 0.666 | 0.668 | 0.628 | 0.484 | 0.402 | 0.436
|  | 10 | 0.546 | 0.680 | 0.654 | 0.648 | 0.646 | 0.558 | 0.488 | 0.504 | 0.336
|  | 11 | 0.548 | 0.706 | 0.682 | 0.664 | 0.714 | 0.594 | 0.576 | 0.578 | 0.450
|  | 12 | 0.634 | 0.684 | 0.704 | 0.678 | 0.640 | 0.534 | 0.608 | 0.664 | 0.472
|  | 1 | 0.688 | 0.722 | 0.642 | 0.696 | 0.684 | 0.698 | 0.684 | 0.678 | 0.646
|  | 2 | 0.656 | 0.718 | 0.684 | 0.686 | 0.642 | 0.642 | 0.606 | 0.626 | 0.564
|  | 3 | 0.670 | 0.692 | 0.674 | 0.692 | 0.652 | 0.634 | 0.632 | 0.680 | 0.610
|  | 4 | 0.696 | 0.656 | 0.654 | 0.724 | 0.664 | 0.682 | 0.592 | 0.658 | 0.610
|  | 5 | 0.730 | 0.708 | 0.652 | 0.732 | 0.680 | 0.642 | 0.630 | 0.664 | 0.652
1 | 50 | 6 | 0.636 | 0.656 | 0.650 | 0.708 | 0.730 | 0.682 | 0.668 | 0.706 | 0.666
|  | 7 | 0.770 | 0.660 | 0.608 | 0.688 | 0.678 | 0.650 | 0.706 | 0.766 | 0.764
|  | 8 | 0.658 | 0.662 | 0.666 | 0.676 | 0.696 | 0.686 | 0.594 | 0.664 | 0.606
|  | 9 | 0.604 | 0.700 | 0.626 | 0.634 | 0.702 | 0.658 | 0.672 | 0.618 | 0.608
|  | 10 | 0.636 | 0.694 | 0.660 | 0.698 | 0.678 | 0.662 | 0.600 | 0.652 | 0.438
|  | 11 | 0.640 | 0.734 | 0.670 | 0.708 | 0.668 | 0.700 | 0.668 | 0.714 | 0.540
|  | 12 | 0.702 | 0.704 | 0.698 | 0.694 | 0.678 | 0.698 | 0.646 | 0.670 | 0.586
|  | 1 | 0.586 | 0.578 | 0.644 | 0.668 | 0.670 | 0.592 | 0.650 | 0.542 | 0.620
|  | 2 | 0.636 | 0.676 | 0.698 | 0.674 | 0.662 | 0.564 | 0.498 | 0.480 | 0.452
|  | 3 | 0.694 | 0.672 | 0.700 | 0.654 | 0.662 | 0.550 | 0.566 | 0.566 | 0.506
|  | 4 | 0.796 | 0.596 | 0.642 | 0.706 | 0.690 | 0.612 | 0.488 | 0.602 | 0.496
|  | 5 | 0.766 | 0.610 | 0.598 | 0.648 | 0.600 | 0.584 | 0.634 | 0.624 | 0.596
7 | 20 | 6 | 0.642 | 0.602 | 0.658 | 0.726 | 0.706 | 0.656 | 0.688 | 0.670 | 0.672
|  | 7 | 0.726 | 0.500 | 0.588 | 0.654 | 0.624 | 0.614 | 0.580 | 0.618 | 0.548
|  | 8 | 0.562 | 0.660 | 0.666 | 0.658 | 0.680 | 0.602 | 0.530 | 0.470 | 0.472
|  | 9 | 0.550 | 0.580 | 0.608 | 0.678 | 0.674 | 0.622 | 0.484 | 0.398 | 0.416
|  | 10 | 0.556 | 0.684 | 0.662 | 0.660 | 0.660 | 0.560 | 0.490 | 0.500 | 0.336
|  | 11 | 0.564 | 0.708 | 0.644 | 0.676 | 0.710 | 0.602 | 0.572 | 0.570 | 0.454
|  | 12 | 0.680 | 0.694 | 0.698 | 0.694 | 0.644 | 0.566 | 0.610 | 0.654 | 0.496
|  | 1 | 0.692 | 0.702 | 0.644 | 0.688 | 0.702 | 0.690 | 0.680 | 0.678 | 0.664
|  | 2 | 0.640 | 0.724 | 0.686 | 0.690 | 0.640 | 0.634 | 0.608 | 0.632 | 0.556
|  | 3 | 0.664 | 0.700 | 0.674 | 0.682 | 0.646 | 0.656 | 0.628 | 0.672 | 0.600
|  | 4 | 0.698 | 0.664 | 0.648 | 0.734 | 0.682 | 0.686 | 0.586 | 0.666 | 0.606
|  | 5 | 0.730 | 0.710 | 0.636 | 0.726 | 0.686 | 0.642 | 0.636 | 0.668 | 0.656
7 | 50 | 6 | 0.642 | 0.652 | 0.664 | 0.706 | 0.742 | 0.700 | 0.676 | 0.716 | 0.702
|  | 7 | 0.774 | 0.622 | 0.610 | 0.694 | 0.700 | 0.666 | 0.704 | 0.758 | 0.744
|  | 8 | 0.658 | 0.662 | 0.658 | 0.678 | 0.702 | 0.680 | 0.600 | 0.656 | 0.610
|  | 9 | 0.618 | 0.690 | 0.622 | 0.638 | 0.704 | 0.654 | 0.690 | 0.640 | 0.612
|  | 10 | 0.634 | 0.682 | 0.642 | 0.688 | 0.678 | 0.642 | 0.600 | 0.654 | 0.436
|  | 11 | 0.646 | 0.730 | 0.668 | 0.706 | 0.668 | 0.718 | 0.670 | 0.718 | 0.552
|  | 12 | 0.700 | 0.708 | 0.714 | 0.700 | 0.678 | 0.720 | 0.644 | 0.664 | 0.572
Table 1: Simulated coverage of 68%-confidence intervals of seasonal structural impulse response of IP ($y^{ip}$), INF $(\pi)$ and FFR ($i$) $k=12$ months after a shock to aggregate supply $(as)$, monetary policy ($mp$) and aggregate demand $(ad)$. Structural schocks are generated using GARCH specification G0. $b$ | $N$ | $s$ | $as\to y^{ip}$ | $as\to\pi$ | $as\to i$ | $mp\to y^{ip}$ | $mp\to\pi$ | $mp\to i$ | $ad\to y^{ip}$ | $ad\to\pi$ | $ad\to i$
---|---|---|---|---|---|---|---|---|---|---|---
|  | 1 | 0.546 | 0.590 | 0.670 | 0.686 | 0.664 | 0.510 | 0.614 | 0.556 | 0.570
|  | 2 | 0.592 | 0.656 | 0.662 | 0.668 | 0.622 | 0.506 | 0.482 | 0.438 | 0.422
|  | 3 | 0.650 | 0.620 | 0.638 | 0.648 | 0.616 | 0.526 | 0.552 | 0.536 | 0.402
|  | 4 | 0.760 | 0.580 | 0.600 | 0.666 | 0.652 | 0.558 | 0.432 | 0.594 | 0.394
|  | 5 | 0.708 | 0.546 | 0.554 | 0.650 | 0.580 | 0.504 | 0.610 | 0.636 | 0.514
1 | 20 | 6 | 0.614 | 0.502 | 0.602 | 0.684 | 0.692 | 0.646 | 0.698 | 0.622 | 0.664
|  | 7 | 0.712 | 0.534 | 0.520 | 0.650 | 0.610 | 0.542 | 0.560 | 0.606 | 0.530
|  | 8 | 0.526 | 0.612 | 0.646 | 0.666 | 0.608 | 0.576 | 0.490 | 0.430 | 0.422
|  | 9 | 0.548 | 0.646 | 0.594 | 0.624 | 0.688 | 0.548 | 0.434 | 0.396 | 0.388
|  | 10 | 0.534 | 0.690 | 0.600 | 0.672 | 0.638 | 0.520 | 0.476 | 0.494 | 0.280
|  | 11 | 0.536 | 0.690 | 0.640 | 0.698 | 0.628 | 0.532 | 0.522 | 0.508 | 0.364
|  | 12 | 0.584 | 0.668 | 0.652 | 0.682 | 0.654 | 0.494 | 0.594 | 0.618 | 0.460
|  | 1 | 0.588 | 0.642 | 0.664 | 0.726 | 0.660 | 0.700 | 0.646 | 0.660 | 0.652
|  | 2 | 0.628 | 0.662 | 0.654 | 0.680 | 0.650 | 0.666 | 0.590 | 0.596 | 0.510
|  | 3 | 0.640 | 0.654 | 0.634 | 0.678 | 0.656 | 0.682 | 0.578 | 0.602 | 0.492
|  | 4 | 0.664 | 0.606 | 0.606 | 0.688 | 0.664 | 0.696 | 0.520 | 0.666 | 0.568
|  | 5 | 0.716 | 0.628 | 0.624 | 0.660 | 0.668 | 0.658 | 0.648 | 0.686 | 0.584
1 | 50 | 6 | 0.610 | 0.618 | 0.666 | 0.706 | 0.668 | 0.662 | 0.640 | 0.712 | 0.656
|  | 7 | 0.728 | 0.594 | 0.558 | 0.656 | 0.654 | 0.624 | 0.610 | 0.702 | 0.670
|  | 8 | 0.636 | 0.634 | 0.614 | 0.656 | 0.644 | 0.656 | 0.540 | 0.578 | 0.514
|  | 9 | 0.550 | 0.684 | 0.598 | 0.686 | 0.666 | 0.710 | 0.614 | 0.580 | 0.556
|  | 10 | 0.606 | 0.680 | 0.620 | 0.712 | 0.664 | 0.648 | 0.594 | 0.636 | 0.376
|  | 11 | 0.600 | 0.706 | 0.670 | 0.716 | 0.652 | 0.662 | 0.638 | 0.662 | 0.504
|  | 12 | 0.568 | 0.666 | 0.616 | 0.672 | 0.676 | 0.694 | 0.622 | 0.652 | 0.508
|  | 1 | 0.540 | 0.592 | 0.654 | 0.690 | 0.670 | 0.534 | 0.622 | 0.566 | 0.566
|  | 2 | 0.600 | 0.666 | 0.646 | 0.686 | 0.632 | 0.496 | 0.494 | 0.436 | 0.412
|  | 3 | 0.666 | 0.620 | 0.660 | 0.674 | 0.622 | 0.520 | 0.564 | 0.542 | 0.414
|  | 4 | 0.780 | 0.598 | 0.612 | 0.686 | 0.658 | 0.584 | 0.446 | 0.584 | 0.384
|  | 5 | 0.726 | 0.548 | 0.558 | 0.664 | 0.584 | 0.510 | 0.632 | 0.646 | 0.548
7 | 20 | 6 | 0.620 | 0.526 | 0.628 | 0.708 | 0.700 | 0.632 | 0.710 | 0.630 | 0.688
|  | 7 | 0.722 | 0.550 | 0.526 | 0.660 | 0.614 | 0.580 | 0.572 | 0.594 | 0.530
|  | 8 | 0.542 | 0.628 | 0.658 | 0.676 | 0.602 | 0.576 | 0.504 | 0.434 | 0.422
|  | 9 | 0.554 | 0.634 | 0.596 | 0.648 | 0.710 | 0.562 | 0.446 | 0.402 | 0.386
|  | 10 | 0.556 | 0.706 | 0.620 | 0.690 | 0.632 | 0.536 | 0.482 | 0.498 | 0.288
|  | 11 | 0.544 | 0.698 | 0.638 | 0.710 | 0.648 | 0.540 | 0.532 | 0.518 | 0.364
|  | 12 | 0.578 | 0.682 | 0.672 | 0.708 | 0.662 | 0.492 | 0.606 | 0.636 | 0.460
|  | 1 | 0.614 | 0.642 | 0.654 | 0.724 | 0.680 | 0.708 | 0.660 | 0.674 | 0.666
|  | 2 | 0.644 | 0.664 | 0.678 | 0.670 | 0.676 | 0.678 | 0.596 | 0.602 | 0.526
|  | 3 | 0.644 | 0.654 | 0.638 | 0.678 | 0.674 | 0.698 | 0.586 | 0.614 | 0.514
|  | 4 | 0.674 | 0.606 | 0.584 | 0.716 | 0.686 | 0.692 | 0.534 | 0.658 | 0.584
|  | 5 | 0.712 | 0.658 | 0.636 | 0.664 | 0.668 | 0.672 | 0.634 | 0.676 | 0.602
7 | 50 | 6 | 0.642 | 0.654 | 0.690 | 0.704 | 0.676 | 0.652 | 0.648 | 0.714 | 0.676
|  | 7 | 0.732 | 0.602 | 0.554 | 0.652 | 0.662 | 0.630 | 0.608 | 0.732 | 0.672
|  | 8 | 0.656 | 0.652 | 0.638 | 0.662 | 0.666 | 0.648 | 0.552 | 0.600 | 0.524
|  | 9 | 0.570 | 0.684 | 0.620 | 0.692 | 0.672 | 0.708 | 0.626 | 0.584 | 0.544
|  | 10 | 0.614 | 0.676 | 0.636 | 0.698 | 0.670 | 0.650 | 0.578 | 0.646 | 0.396
|  | 11 | 0.608 | 0.696 | 0.692 | 0.708 | 0.660 | 0.672 | 0.644 | 0.666 | 0.506
|  | 12 | 0.570 | 0.682 | 0.650 | 0.684 | 0.696 | 0.698 | 0.632 | 0.660 | 0.512
Table 2: Simulated coverage of 68%-confidence intervals of seasonal structural
impulse response of industrial production ($y^{ip}$), inflation $(\pi)$ and
federal funds rate ($i$) $k=12$ months after a shock to aggregate supply
$(as)$, monetary policy ($mp$) and aggregate demand $(ad)$. Structural schocks
are generated using GARCH specification G3.

We simulate a total of $M=500$ time series of length $N=20,50,100$ resulting
in $T=240,600,1200$. In each Monte Carlo iteration, we ensure that the
simulated FFR does not fall below zero by setting negative values to its
natural lower bound 0. Further, we bootstrap the structural impulse responses
by using the residual-based seasonal block bootstrap for PVAR processes
discussed in Section 5.1.1 with block sizes $b\in\\{1,7\\},$ in order to
obtain approximations of the confidence intervals. Note that for $b=1$, the
residual-based seasonal block bootstrap collapses to a seasonal independent
variant. The confidence intervals are constructed according to (14). The
nominal coverage rate is 68% and we use $B=500$ bootstrap repetitions.

Tables 1 and 2 show the simulated coverage rates for lag $k=12$ and
$N\in\\{20,50\\}$ for GARCH specifications G0 and G3, where the seasonal block
bootstrap ($b=7$) and its independent version ($b=1$) are used for each
specification. The tables for lag $k=0$, $N=100$ and for the other GARCH
specifications G1 and G2 can be found in supplementary appendix F,
respectively. Due to a simplified presentation, we only provide coverage rates
for lags $k\in\\{0,12\\}$. For other lags, the coverage rate is not
particularly different. From Tables 1 and 2, it can be clearly seen that the
simulated coverage rates approach the nominal coverage rate of 68% for
increasing $N$. When we consider Table 1, very low efficiency losses are
observed for $N\in\\{50,100\\}$ when using the seasonal block bootstrap
$(b=7)$. For $N=20$, the seasonal block bootstrap $(b=7)$ achieves the
coverage slightly better. This loss of efficiency is due to the fact that the
seasonal block bootstrap is not required for the i.i.d. case G0. When we
consider Table 2, the simulated coverage rates approach the nominal coverage
rate somewhat more closely when using the seasonal block bootstrap. We also
find this result for cases G1 and G2. This is because the seasonal independent
bootstrap is not able to capture the non-linear dependencies of the GARCH
errors. While the simulated coverage rate for most impulse responses is very
close to the correct coverage rate, we see only minor deviations from the
nominal coverage in the impulse responses of the FFR after an aggregate demand
shock. Compared to the standard percentile, Hall’s percentile or Hall’s
studentized interval, the method for the construction of confidence intervals
proposed in (14) achieves the correct coverage slightly better, especially for
the impulse response of the FFR after an aggregate demand shock.

##  7 Conclusion

When dealing with seasonal data, deseasonalization using standard seasonal
adjustment procedures is the standard methodology nowadays, even though
seasonal adjustment procedures may distort the structure of the data and
remove useful information. In this paper, we offer an alternative approach to
analyze not seasonally adjusted macroeconomic data by structural models.
Instead of first seasonal adjusting the data and then fitting a structural
VAR, we propose to fit a structural periodic VAR directly on not seasonally
adjusted macroeconomic raw data. Moreover, we provide a PVAR representation
that enables linearly restricted estimation of PVAR models under more general
linear constraints than Ursu & Duchesne (2009), which allows to impose also
constraints _across seasons_. Further, we show consistency and asymptotic
normality of (constrained) least squares estimators for the PVAR coefficients
and structural impulse responses under periodic weak white noise assumption.
For standard error estimation and confidence interval construction, we also
introduce residual-based bootstrap methods for PVARs and prove bootstrap
consistency.

The empirical results of Section 6 illustrate the analysis of not seasonally
adjusted industrial production, CPI inflation and federal funds rate using a
linearly constrained structural PVAR. We identify three different shocks, i.e.
monetary policy shock, aggregate supply and demand shock based on a
combination of short- and long-run restrictions following Peersman (2005). We
determine seasonal structural impulse responses based on the linearly
restricted structural PVAR and compare the results with those from a standard
SVAR analysis on the same US macro variables in seasonally adjusted form. We
find that monetary policy shocks do not seem to have a periodic effect on the
macroeconomic variables. However, we find some significant periodic effects
caused by aggregate supply and demand shocks, especially in July, which may
indicate a cetain holiday effect.

Moreover, we do not find pronounced efficiency losses of suitably constrained
structural PVARs compared to seasonally adjusted SVARs, but even efficiency
gains in individual seasons. In total, we clearly see that useful insights
into the dynamics of the macro variables are lost, if we seasonally adjust the
data. In the simulation study, we showed that the simulated coverage rates
constructed by the proposed residual-based bootstrap methods for PVARs achieve
the nominal coverage rate.

## References

  * (1)
  * Aknouche (2007) Aknouche, A. (2007), ‘Causality conditions and autocovariance calculations in PVAR models’, Journal of Statistical Computation and Simulation 77(9), 769–780.
  * Aknouche (2017) Aknouche, A. (2017), ‘Periodic autoregressive stochastic volatility’, Statistical Inference for Stochastic Processes 20, 139–177.
  * Aknouche & Bibi (2009) Aknouche, A. & Bibi, A. (2009), ‘Quasi-maximum likelihood estimation of periodic GARCH and periodic ARMA-GARCH processes’, Journal of Time Series Analysis 30(1), 19–46.
  * Bell (2012) Bell, W. R. (2012), ‘Unit root properties of seasonal adjustment and related filters’, Journal of Official Statistics 28(3), 441–461.
  * Bell & Hillmer (1984) Bell, W. R. & Hillmer, S. C. (1984), ‘Issues involved with the seasonal adjustment of economic time series’, Journal of Business & Economic Statistics 2(4), 291–320.
  * Bernanke & Blinder (1992) Bernanke, B. & Blinder, A. S. (1992), ‘The federal funds rate and the transmission of monetary policy’, American Economic Review 82(4), 901–921.
  * Bernanke et al. (1997) Bernanke, B. S., Gertler, M., Watson, M., Sims, C. A. & Friedman, B. M. (1997), ‘Systematic monetary policy and the effects of oil price shocks’, Brookings papers on economic activity 1997(1), 91–157.
  * Bertail & Dudek (2022) Bertail, P. & Dudek, A. E. (2022), ‘Optimal choice of bootstrap block length for periodically correlated time series’, hal-03472199v2 .
  * Bibi & Ghezal (2016) Bibi, A. & Ghezal, A. (2016), ‘On periodic time-varying bilinear processes: Structure and asymptotic inference’, Statistical Methods & Applications 25, 395–420.
  * Blanchard & Quah (1989) Blanchard, O. J. & Quah, D. (1989), ‘The dynamic effects of aggregate demand and supply disturbances’, 79(4), 655–673.
  * Bollerslev & Ghysels (1996) Bollerslev, T. & Ghysels, E. (1996), ‘Periodic autoregressive conditional heteroscedasticity’, Journal of Business & Economic Statistics 14(2), 139–151.
  * Boubacar Maïnassara & Ursu (2023) Boubacar Maïnassara, Y. & Ursu, E. (2023), ‘Estimating weak periodic vector autoregressive time series’, TEST 32, 958–997.
  * Brockwell & Davies (1991) Brockwell, P. & Davies, R. (1991), Time Series: Theory and Methods, Springer Verlag.
  * Brüggemann et al. (2016) Brüggemann, R., Jentsch, C. & Trenkler, C. (2016), ‘Inference in VARs with conditional heteroskedasticity of unknown form’, Journal of Econometrics 191(1), 69–85.
  * Canova & Hansen (1995) Canova, F. & Hansen, B. E. (1995), ‘Are seasonal patterns constant over time? A test for seasonal stability’, Journal of Business & Economic Statistics 13(3), 237–252.
  * Consolvo & Lunsford (2019) Consolvo, V. & Lunsford, K. G. (2019), ‘Residual seasonality in GDP growth remains after latest BEA improvements’, Economic Commentary (2019-05).
  * Del Barrio Castro & Osborn (2004) Del Barrio Castro, T. & Osborn, D. R. (2004), ‘The consequences of seasonal adjustment for periodic autoregressive processes’, The Econometrics Journal 7(2), 307–321.
  * Doppelt (2021) Doppelt, R. (2021), ‘Should macroeconomists use seasonally adjusted time series? Structural identification and bayesian estimation in seasonal vector autoregressions’, Manuscript, Michigan State University. https://www2. gwu. edu/~ forcpgm/Doppelt_Seasonality_2021. pdf .
  * Dudek et al. (2014) Dudek, A. E., Leśkow, J., Paparoditis, E. & Politis, D. N. (2014), ‘A generalized block bootstrap for seasonal time series’, Journal of Time Series Analysis 35(2), 89–114.
  * Dudek et al. (2016) Dudek, A. E., Paparoditis, E. & Politis, D. N. (2016), ‘Generalized seasonal tapered block bootstrap’, Statistics & Probability Letters 115, 27–35.
  * Findley et al. (2017) Findley, D. F., Lytras, D. P. & McElroy, T. S. (2017), ‘Detecting seasonality in seasonally adjusted monthly time series’.  Census Bureau Research Report: RRS2017/03.
  * Findley et al. (1998) Findley, D. F., Monsell, B. C., Bell, W. R., Otto, M. C. & Chen, B.-C. (1998), ‘New capabilities and methods of the X-12-ARIMA seasonal-adjustment program’, Journal of Business & Economic Statistics 16(2), 127–152.
  * Francq & Raïssi (2007) Francq, C. & Raïssi, H. (2007), ‘Multivariate portmanteau test for autoregressive models with uncorrelated but nonindependent errors’, Journal of Time Series Analysis 28(3), 454–470.
  * Francq et al. (2011) Francq, C., Roy, R. & Saidi, A. (2011), ‘Asymptotic properties of weighted least squares estimation in weak PARMA models’, Journal of Time Series Analysis 32(6), 699–723.
  * Franses (1996) Franses, P. H. (1996), ‘Recent advances in modelling seasonality’, Journal of Economic Surveys 10(3), 299–345.
  * Franses & Paap (2004) Franses, P. H. & Paap, R. (2004), Periodic time series models, Oxford University Press.
  * Franses & Van Dijk (2005) Franses, P. H. & Van Dijk, D. (2005), ‘The forecasting performance of various models for seasonality and nonlinearity for quarterly industrial production’, International Journal of Forecasting 21(1), 87–102.
  * Gali (1992) Gali, J. (1992), ‘How well does the IS-LM model fit postwar US data?’, The Quarterly Journal of Economics 107(2), 709–738.
  * Gerlach & Smets (1995) Gerlach, S. & Smets, F. (1995), ‘The monetary transmission mechanism: Evidence from the G-7 countries’, Centre for Economic Research Discussion paper No. 1219 .
  * Gersovitz & MacKinnon (1978) Gersovitz, M. & MacKinnon, J. G. (1978), ‘Seasonality in regression: An application of smoothness priors’, Journal of the American Statistical Association 73(362), 264–273.
  * Ghysels et al. (1996) Ghysels, E., Granger, C. W. & Siklos, P. L. (1996), ‘Is seasonal adjustment a linear or nonlinear data-filtering process?’, Journal of Business & Economic Statistics 14(3), 374–386.
  * Ghysels & Osborn (2001) Ghysels, E. & Osborn, D. R. (2001), The econometric analysis of seasonal time series, Cambridge University Press.
  * Ghysels et al. (2006) Ghysels, E., Osborn, D. R. & Rodrigues, P. M. (2006), ‘Forecasting seasonal time series’, Handbook of Economic Forecasting 1, 659–711.
  * Gladyshev (1961) Gladyshev, E. (1961), ‘Periodically correlated random sequence’, Soviet. Math. 2, 385–388.
  * Gómez & Maravall Herrero (1996) Gómez, V. & Maravall Herrero, A. (1996), Programs TRAMO and SEATS: Instructions for the user (beta version: September 1996), Banco de España. Servicio de Estudios.
  * Gonçalves & Kilian (2007) Gonçalves, S. & Kilian, L. (2007), ‘Asymptotic and bootstrap inference for AR($\infty$) processes with conditional heteroskedasticity’, Econometric Reviews 26(6), 609–641.
  * Hylleberg (1992) Hylleberg, S. (1992), Modelling seasonality, Oxford University Press.
  * Hylleberg et al. (1993) Hylleberg, S., Jørgensen, C. & Sørensen, N. K. (1993), ‘Seasonality in macroeconomic time series’, Empirical Economics 18, 321–335.
  * Jentsch & Lunsford (2019) Jentsch, C. & Lunsford, K. G. (2019), ‘The dynamic effects of personal and corporate income tax changes in the United States: Comment’, American Economic Review 109(7), 2655–2678.
  * Jentsch & Lunsford (2022) Jentsch, C. & Lunsford, K. G. (2022), ‘Asymptotically valid bootstrap inference for proxy SVARs’, Journal of Business & Economic Statistics 40(4), 1876–1891.
  * Jones & Brelsford (1967) Jones, R. H. & Brelsford, W. M. (1967), ‘Time series with periodic structure’, Biometrika 54(3-4), 403–408.
  * Kilian (2009) Kilian, L. (2009), ‘Not all oil price shocks are alike: Disentangling demand and supply shocks in the crude oil market’, American Economic Review 99(3), 1053–1069.
  * Kilian & Lütkepohl (2017) Kilian, L. & Lütkepohl, H. (2017), Structural vector autoregressive analysis, Cambridge University Press.
  * Kilian & Murphy (2012) Kilian, L. & Murphy, D. P. (2012), ‘Why agnostic sign restrictions are not enough: Understanding the dynamics of oil market VAR models’, Journal of the European Economic Association 10(5), 1166–1188.
  * Lund & Basawa (2000) Lund, R. & Basawa, I. (2000), ‘Recursive prediction and likelihood evaluation for periodic ARMA models’, Journal of Time Series Analysis 21(1), 75–93.
  * Lunsford (2017) Lunsford, K. G. (2017), ‘Lingering residual seasonality in GDP growth’, Economic Commentary (2017-06).
  * Lütkepohl (2005) Lütkepohl, H. (2005), New introduction to multiple time series analysis, Springer Science & Business Media.
  * Monsell (2007) Monsell, B. C. (2007), The X-13A-S seasonal adjustment program, in ‘Proceedings of the 2007 Federal Committee On Statistical Methodology Research Conference. URL http://www. fcsm. gov/07papers/Monsell. II-B. pdf’, p. 515.
  * Nordman (2009) Nordman, D. J. (2009), ‘A note on the stationary bootstrap’s variance’, The Annals of Statistics 37, 359–370.
  * Novales & de Fruto (1997) Novales, A. & de Fruto, R. F. (1997), ‘Forecasting with periodic models a comparison with time invariant coefficient models’, International Journal of Forecasting 13(3), 393–405.
  * Osborn (1988) Osborn, D. R. (1988), ‘Seasonality and habit persistence in a life cycle model of consumption’, Journal of Applied Econometrics 3, 255–266.
  * Osborn & Smith (1989) Osborn, D. R. & Smith, J. P. (1989), ‘The performance of periodic autoregressive models in forecasting seasonal UK consumption’, Journal of Business & Economic Statistics 7(1), 117–127.
  * Pagano (1978) Pagano, M. (1978), ‘On periodic and multiple autoregressions’, The Annals of Statistics 6, 1310–1317.
  * Peersman (2005) Peersman, G. (2005), ‘What caused the early millennium slowdown? Evidence based on vector autoregressions’, Journal of Applied Econometrics 20(2), 185–207.
  * Rudebusch et al. (2015) Rudebusch, G. D., Wilson, D. & Mahedy, T. (2015), ‘The puzzle of weak first-quarter GDP growth’, Puzzle 2015, 16.
  * Shiskin et al. (1967) Shiskin, J., Young, A. & Musgrave, J. (1967), ‘The X-11 variant of the Census II method seasonal adjustment program’.  Technical Paper No. 15, US Department of Commerce, Bureau of Economic Analysis.
  * Troutman (1979) Troutman, B. M. (1979), ‘Some results in periodic autoregression’, Biometrika 66(2), 219–228.
  * Ula (1990) Ula, T. A. (1990), ‘Periodic covariance stationarity of multivariate periodic autoregressive moving average processes’, Water Resources Research 26(5), 855–861.
  * Ula (1993) Ula, T. A. (1993), ‘Forecasting of multivariate periodic autoregressive moving-average processes’, Journal of Time Series Analysis 14(6), 645–657.
  * Ursu & Duchesne (2009) Ursu, E. & Duchesne, P. (2009), ‘On modelling and diagnostic checking of vector periodic autoregressive time series models’, Journal of Time Series Analysis 30(1), 70–96.

Online Supplemental Appendices to

Structural Periodic Vector Autoregressions

Daniel Dzikowski§§§TU Dortmund University, Department of Statistics, D-44221
Dortmund, Germany<EMAIL_ADDRESS>corresponding author
TU Dortmund University

Carsten Jentsch¶¶¶TU Dortmund University, Department of Statistics, D-44221
Dortmund, Germany<EMAIL_ADDRESS>
TU Dortmund University

August 28, 2024

##  Appendix A Proof of Theorem 3.1

In order to prove the joint CLT of the PVAR estimators under Assumption 3.1,
some pre-work has to be done. At first, for $a,b,c\in\mathbb{Z}$, we define
$\kappa_{a,b}(s_{1},s_{2})$, $\tau_{a,b,c}(s_{1},s_{2})$,
$s_{1},s_{2}=1,\dots,S$ and $\tau_{a,b,c}$ as in (11) and (12). As given in
(10), the multivariate LS estimator under general linear restrictions
$\beta=R\gamma+r$ is given by:

$\displaystyle\widehat{\gamma}=\left[R^{\prime}\\{XX^{\prime}\otimes\bm{I}_{m}\\}R\right]^{-1}R^{\prime}\\{X\otimes\bm{I}_{m}\\}\mathbb{z}.$

Plugging in $\mathbb{z}=\\{X^{\prime}\otimes\bm{I}_{m}\\}R\gamma+e$ gives

$\displaystyle\widehat{\gamma}-\gamma=\left[R^{\prime}\\{XX^{\prime}\otimes\bm{I}_{m}\\}R\right]^{-1}R^{\prime}\\{X\otimes\bm{I}_{m}\\}e=\left[R^{\prime}\\{XX^{\prime}\otimes\bm{I}_{m}\\}R\right]^{-1}R^{\prime}\;vec\\{EX^{\prime}\\}$

where

$vec\\{EX^{\prime}\\}=\sum\limits_{n=0}^{N-1}\begin{pmatrix}\epsilon_{Sn+1}\
\\\ vec\\{\epsilon_{Sn+1}y_{Sn}^{\prime}\\}\\\ \vdots\\\
vec\\{\epsilon_{Sn+1}y_{Sn-p(1)+1}^{\prime}\\}\\\ \epsilon_{Sn+2}\\\
vec\\{\epsilon_{Sn+2}y_{Sn+1}^{\prime}\\}\\\ \vdots\\\
vec\\{\epsilon_{Sn+2}y_{Sn-p(2)+2}^{\prime}\\}\\\ \vdots\\\
\epsilon_{Sn+S-1}\\\ vec\\{\epsilon_{Sn+S}y_{Sn+S-1}^{\prime}\\}\\\ \vdots\\\
vec\\{\epsilon_{Sn+S}y_{Sn+S-p(S)}^{\prime}\\}\end{pmatrix}.$

Due to Assumption 3.1 $(i)$, the moving-average representation
$y_{Sn+s}=\mu(s)+\sum\limits_{k=0}^{\infty}\Phi_{k}(s)\epsilon_{Sn+s-k}$ of
the PVAR($\bm{p}$) process exists and can be used to represent $X_{n}(s)$ as

$\displaystyle X_{n}(s)$ $\displaystyle=\begin{pmatrix}1\\\ y_{Sn+s-1}\\\
\vdots\\\ y_{Sn+s-p(s)}\end{pmatrix}=\begin{pmatrix}1\\\
\mu(s-1)+\sum\limits_{k=0}^{\infty}\Phi_{k}(s-1)\epsilon_{Sn+s-1-k}\\\
\vdots\\\
\mu(s-p(s))+\sum\limits_{k=0}^{\infty}\Phi_{k}(s-p(s))\epsilon_{Sn+s-p(s)-k}\end{pmatrix}$
$\displaystyle=\begin{pmatrix}1\\\
\sum\limits_{k=1}^{\infty}\begin{pmatrix}\Phi_{k-1}(s-1)[\epsilon_{Sn+s-k}+\nu(s-k)]\\\
\vdots\\\
\Phi_{k-1}(s-p(s))[\epsilon_{Sn+s-p(s)+1-k}+\nu(s-p(s)+1-k)]\end{pmatrix}\end{pmatrix}$
$\displaystyle=\begin{pmatrix}1\\\
\sum\limits_{k=1}^{\infty}\begin{pmatrix}\Phi_{k-1}(s-1)[\epsilon_{Sn+s-k}+\nu(s-k)]\\\
\vdots\\\
\Phi_{k-p(s)}(s-p(s))[\epsilon_{Sn+s-k}+\nu(s-k)]\end{pmatrix}\end{pmatrix}$
$\displaystyle=\begin{pmatrix}1\\\
\sum\limits_{k=1}^{\infty}C_{k}(s)[\epsilon_{Sn+s-k}+\nu(s-k)]\end{pmatrix},$

where

$C_{k}(s)=\begin{pmatrix}\Phi_{k-1}(s-1)\\\ \vdots\\\
\Phi_{k-p(s)}(s-p(s))\end{pmatrix}\in\mathbb{R}^{mp(s)\times m}.$

Note that $\Phi_{k-p(s)}(s)=\bm{0}$ for $k<p(s),s=1,\dots,S$. Hence, for
$vec\\{EX^{\prime}\\}$, we get

$\displaystyle vec\\{EX^{\prime}\\}$
$\displaystyle=\begin{pmatrix}\bm{I}_{m}&\bm{0}_{m\times
m^{2}}&\dots&\dots&\bm{0}\\\ \bm{0}_{m^{2}p(1)\times
m}&\sum\limits_{k=1}^{\infty}(C_{k}(1)\otimes\bm{I}_{m})&\ddots&&\vdots\\\
\vdots&\ddots&\ddots&\ddots&\vdots\\\
\vdots&&\ddots&\bm{I}_{m}&\bm{0}_{m\times m^{2}}\\\
\bm{0}&\dots&\dots&\bm{0}_{m^{2}p(S)\times
m}&\sum\limits_{k=1}^{\infty}(C_{k}(S)\otimes\bm{I}_{m})\end{pmatrix}$
$\displaystyle\times\sum\limits_{n=0}^{N-1}\begin{pmatrix}\epsilon_{Sn+1}\\\
vec\\{\epsilon_{Sn+1}[\epsilon_{Sn+1-k}^{\prime}+\nu(1-k)^{\prime}]\\}\\\
\epsilon_{Sn+2}\\\
vec\\{\epsilon_{Sn+2}[\epsilon_{Sn+2-k}^{\prime}+\nu(2-k)^{\prime}]\\}\\\
\vdots\\\ \epsilon_{Sn+S}\\\
vec\\{\epsilon_{Sn+S}[\epsilon_{Sn+S-k}^{\prime}+\nu(S-k)^{\prime}]\\}\end{pmatrix}.$

Further, using
$\sqrt{N}\bigl{(}\widehat{\beta}_{res}-\beta\bigr{)}=\sqrt{N}R\bigl{(}\widehat{\gamma}-\gamma\bigr{)}$,
we get the deviation form of $\beta$, that is

$\displaystyle\sqrt{N}\bigl{(}\widehat{\beta}_{res}-\beta\bigr{)}$
$\displaystyle=R\left[R^{\prime}\\{\frac{1}{N}XX^{\prime}\otimes\bm{I}_{m}\\}R\right]^{-1}R^{\prime}$
$\displaystyle\times\begin{pmatrix}\bm{I}_{m}&\bm{0}_{m\times
m^{2}}&\dots&\dots&\bm{0}\\\ \bm{0}_{m^{2}p(1)\times
m}&\sum\limits_{k=1}^{\infty}(C_{k}(1)\otimes\bm{I}_{m})&\ddots&&\vdots\\\
\vdots&\ddots&\ddots&\ddots&\vdots\\\
\vdots&&\ddots&\bm{I}_{m}&\bm{0}_{m\times m^{2}}\\\
\bm{0}&\dots&\dots&\bm{0}_{m^{2}p(S)\times
m}&\sum\limits_{k=1}^{\infty}(C_{k}(S)\otimes\bm{I}_{m})\end{pmatrix}$ (A.1)
$\displaystyle\times\frac{1}{\sqrt{N}}\sum\limits_{n=0}^{N-1}\begin{pmatrix}\epsilon_{Sn+1}\\\
vec\\{\epsilon_{Sn+1}[\epsilon_{Sn+1-k}^{\prime}+\nu(1-k)^{\prime}]\\}\\\
\epsilon_{Sn+2}\\\
vec\\{\epsilon_{Sn+2}[\epsilon_{Sn+2-k}^{\prime}+\nu(2-k)^{\prime}]\\}\\\
\vdots\\\ \epsilon_{Sn+S}\\\
vec\\{\epsilon_{Sn+S}[\epsilon_{Sn+S-k}^{\prime}+\nu(S-k)^{\prime}]\\}\end{pmatrix}.$

Next, we derive asymptotic properties of the estimator of the periodic
covariance matrix $\Sigma_{\epsilon}(s)$. For this purpose, we define
$\sigma(s)=vech\\{\Sigma_{\epsilon}(s)\\}$,
$\widehat{\sigma}(s)=vech\\{\widehat{\Sigma}_{\epsilon}(s)\\}$ and
$\widetilde{\sigma}(s)=vech\\{\widetilde{\Sigma}_{\epsilon}(s)\\}$. While
$\widehat{\Sigma}_{\epsilon}(s)=\frac{1}{N-k(s)}\sum_{n=0}^{N-1}\widehat{\epsilon}_{Sn+s}\widehat{\epsilon}_{Sn+s}^{\prime}$
and
$\widetilde{\Sigma}_{\epsilon}(s)=\frac{1}{N}\sum_{n=0}^{N-1}\epsilon_{Sn+s}\epsilon_{Sn+s}^{\prime}$.
$\widehat{\Sigma}_{\epsilon}(s)$ is the natural estimator of the seasonal
covariance matrix $\Sigma_{\epsilon}(s)$, $s=1,\dots,S$ based on the PVAR
residuals, $\widetilde{\Sigma}_{\epsilon}(s)$ is the corresponding
(infeasible) estimator based on the non-observable PVAR error terms. The
$vech$-operator stacks the entries on the lower-triangular part of a square
matrix columnwise below each other. Further, let
$\sigma=(\sigma(1)^{\prime},\dots,\sigma(S)^{\prime})^{\prime}$,
$\widehat{\sigma}=(\widehat{\sigma}(1)^{\prime},\dots,\widehat{\sigma}(S)^{\prime})^{\prime}$
and
$\widetilde{\sigma}=(\widetilde{\sigma}(1)^{\prime},\dots,\widetilde{\sigma}(S)^{\prime})^{\prime}$.
By standard arguments, it can be shown easily that

$\sqrt{N}\bigl{(}\widetilde{\sigma}-\widehat{\sigma}\bigr{)}=o_{P}(1).$

Hence, for further calculations, $\widetilde{\sigma}$ can be used and its
deviation form can be expressed by

$\sqrt{N}\bigl{(}\widetilde{\sigma}-\sigma\bigr{)}=\frac{1}{\sqrt{N}}\sum\limits_{n=0}^{N-1}L_{mS}\begin{pmatrix}vec\\{\epsilon_{Sn+1}\epsilon_{Sn+1}^{\prime}\\}-vec\\{\Sigma_{\epsilon}(1)\\}\\\
\vdots\\\
vec\\{\epsilon_{Sn+S}\epsilon_{Sn+S}^{\prime}\\}-vec\\{\Sigma_{\epsilon}(S)\\}\end{pmatrix}.$

Together with (A), the joint deviation form can be represented as

$\displaystyle\sqrt{N}\begin{pmatrix}\widehat{\beta}_{res}-\beta\\\
\widetilde{\sigma}-\sigma\end{pmatrix}$
$\displaystyle=\begin{pmatrix}\widehat{Q}^{\beta}R^{\beta}\times\frac{1}{\sqrt{N}}\sum_{n=0}^{N-1}\begin{pmatrix}\epsilon_{Sn+1}\\\
vec\\{\epsilon_{Sn+1}[\epsilon_{Sn+1-k}^{\prime}+\nu(1-k)^{\prime}]\\}\\\
\epsilon_{Sn+2}\\\
vec\\{\epsilon_{Sn+2}[\epsilon_{Sn+2-k}^{\prime}+\nu(2-k)^{\prime}]\\}\\\
\vdots\\\ \epsilon_{Sn+S}\\\
vec\\{\epsilon_{Sn+S}[\epsilon_{Sn+S-k}^{\prime}+\nu(S-k)^{\prime}]\\}\end{pmatrix}\\\
\\\
\frac{1}{\sqrt{N}}\sum_{n=0}^{N-1}L_{mS}\begin{pmatrix}vec\\{\epsilon_{Sn+1}\epsilon_{Sn+1}^{\prime}\\}-vec\\{\Sigma_{\epsilon}(1)\\}\\\
\vdots\\\
vec\\{\epsilon_{Sn+S}\epsilon_{Sn+S}^{\prime}\\}-vec\\{\Sigma_{\epsilon}(S)\\}\end{pmatrix}\end{pmatrix}$
(A.2) $\displaystyle=A_{q}+(A-A_{q}),$

with

$\displaystyle\widehat{Q}^{\beta}$
$\displaystyle=R\left[R^{\prime}\\{\frac{1}{N}XX^{\prime}\otimes\bm{I}_{m}\\}R\right]^{-1}R^{\prime}$
$\displaystyle R^{\beta}$
$\displaystyle=\begin{pmatrix}\bm{I}_{m}&\bm{0}_{m\times
m^{2}}&\dots&\dots&\bm{0}\\\ \bm{0}_{m^{2}p(1)\times
m}&\sum_{k=1}^{\infty}(C_{k}(1)\otimes\bm{I}_{m})&\ddots&&\vdots\\\
\vdots&\ddots&\ddots&\ddots&\vdots\\\
\vdots&&\ddots&\bm{I}_{m}&\bm{0}_{m\times m^{2}}\\\
\bm{0}&\dots&\dots&\bm{0}_{m^{2}p(S)\times
m}&\sum_{k=1}^{\infty}(C_{k}(S)\otimes\bm{I}_{m})\end{pmatrix},$

where $A$ denotes the right-hand side of (A.2) and $A_{q}$ is defined as the
same, but with $\sum_{k=1}^{\infty}$ replaced by $\sum_{k=1}^{q}$ for some
$q\in\mathbb{N}$. For proving the Central Limit Theorem of the joint deviation
form, we make use of Proposition 6.3.9 in Brockwell & Davies (1991) and it
suffices to show

  * (a)

$A_{q}\overset{d}{\to}\mathcal{N}(0,V_{q})$ as $N\to\infty$

  * (b)

$V_{q}\to V$ as $q\to\infty$

  * (c)

$\forall\delta>0:\;\;\underset{q\to\infty}{\lim}\underset{N\to\infty}{\limsup}P(|A-A_{q}|_{1}>\delta)=0$

In order to prove (a), setting $\widetilde{m}=\frac{m(m+1)}{2}$, we bring
$A_{q}$ to the following form

$\displaystyle
A_{q}=\widehat{Q}_{N}R_{q}\frac{1}{\sqrt{N}}\sum\limits_{n=0}^{N-1}\begin{pmatrix}W_{n,q}^{(1)}\\\
W_{n,q}^{(2)}\end{pmatrix},$

where

$\displaystyle\widehat{Q}_{N}=\begin{pmatrix}R\left[R^{\prime}\\{\frac{1}{N}XX^{\prime}\otimes\bm{I}_{m}\\}R\right]^{-1}R^{\prime}\;\;\;&\bm{0}_{m\sum_{s=1}^{S}(mp(s)+1)\times
S\widetilde{m}}\\\ \bm{0}_{S\widetilde{m}\times
m\sum_{s=1}^{S}(mp(s)+1)}&\bm{I}_{S\widetilde{m}}\end{pmatrix},$
$\displaystyle
R_{q}=\begin{pmatrix}\bm{I}_{m}&\bm{0}&\dots&\bm{0}&\dots&\dots&\dots&\dots&\dots&\bm{0}\\\
\bm{0}&C_{1}(1)\otimes\bm{I}_{m}&\dots&C_{q}(1)\otimes\bm{I}_{m}&\ddots&&&&&\vdots\\\
\bm{0}&\bm{0}&\dots&\bm{0}&\ddots&\ddots&&&&\vdots\\\
\vdots&&&&\ddots&\bm{I}_{m}&\bm{0}&\dots&\bm{0}&\bm{0}\\\
\vdots&&&&&\bm{0}&C_{1}(S)\otimes\bm{I}_{m}&\dots&C_{q}(S)\otimes\bm{I}_{m}&\bm{0}\\\
\bm{0}&\dots&\dots&\dots&\dots&\bm{0}&\bm{0}&\dots&\bm{0}&\bm{I}_{S\widetilde{m}}\end{pmatrix}$

and

$\displaystyle\begin{pmatrix}W_{n,q}^{(1)}\\\
W_{n,q}^{(2)}\end{pmatrix}=\begin{pmatrix}\begin{pmatrix}\epsilon_{Sn+1}\\\
vec\\{\epsilon_{Sn+1}[\epsilon_{Sn}^{\prime}+\nu(0)^{\prime}]\\}\\\ \vdots\\\
vec\\{\epsilon_{Sn+1}[\epsilon_{Sn+1-q}^{\prime}+\nu(1-q)^{\prime}]\\}\\\
\epsilon_{Sn+2}\\\
vec\\{\epsilon_{Sn+2}[\epsilon_{Sn+1}^{\prime}+\nu(1)^{\prime}]\\}\\\
\vdots\\\
vec\\{\epsilon_{Sn+2}[\epsilon_{Sn+2-q}^{\prime}+\nu(2-q)^{\prime}]\\}\\\
\vdots\\\ \epsilon_{Sn+S}\\\
vec\\{\epsilon_{Sn+S}[\epsilon_{Sn+S-1}^{\prime}+\nu(S-1)^{\prime}]\\}\\\
\vdots\\\
vec\\{\epsilon_{Sn+S}[\epsilon_{Sn+S-q}^{\prime}+\nu(S-q)^{\prime}]\\}\end{pmatrix}\\\
L_{mS}\begin{pmatrix}vec\\{\epsilon_{Sn+1}\epsilon_{Sn+1}^{\prime}\\}-vec\\{\Sigma_{\epsilon}(1)\\}\\\
\vdots\\\
vec\\{\epsilon_{Sn+S}\epsilon_{Sn+S}^{\prime}\\}-vec\\{\Sigma_{\epsilon}(S)\\}\end{pmatrix}\end{pmatrix}$

are of dimension
$(m\sum_{s=1}^{S}(mp(s)+1)+S\widetilde{m})\times(m\sum_{s=1}^{S}(mp(s)+1)+S\widetilde{m}),(m\sum_{s=1}^{S}(mp(s)+1)+S\widetilde{m})\times(S(m+qm^{2})+S\widetilde{m})$
and $S(m+qm^{2})+S\widetilde{m}$. Following Brüggemann et al. (2016),
$\frac{1}{N}XX^{\prime}\overset{p}{\to}\Gamma$ as $N\to\infty$ leads to
$\widehat{Q}_{N}\overset{p}{\to}Q$ as $N\to\infty$ , where

$Q=\begin{pmatrix}R\left[R^{\prime}\\{\Gamma\otimes\bm{I}_{m}\\}R\right]^{-1}R^{\prime}\;\;\;&\bm{0}_{m\sum_{s=1}^{S}(mp(s)+1)\times
S\widetilde{m}}\\\ \bm{0}_{S\widetilde{m}\times
m\sum_{s=1}^{S}(mp(s)+1)}&\bm{I}_{S\widetilde{m}}\end{pmatrix}.$

The matrix $\Gamma=E(\frac{1}{N}XX^{\prime})$ is of dimension
$\sum_{s=1}^{S}(mp(s)+1)\times\sum_{s=1}^{S}(mp(s)+1)$ and can be represented
as

$\Gamma=\begin{pmatrix}\Gamma(1)&0&\dots&0\\\ 0&\Gamma(2)&\ddots&\vdots\\\
\vdots&\ddots&\ddots&0\\\ 0&\dots&0&\Gamma(S)\end{pmatrix},$

where $\Gamma(s)$, $s=1,\dots,S$ is given by

$\Gamma(s)=\begin{pmatrix}1&\widetilde{\mu}(s)^{\prime}\\\
\widetilde{\mu}(s)&\widetilde{\mu}(s)\widetilde{\mu}(s)^{\prime}+\sum_{k=1}^{\infty}C_{k}(s)\Sigma_{\epsilon}(s-k)C_{k}^{\prime}(s)\end{pmatrix}\in\mathbb{R}^{(mp(s)+1)\times(mp(s)+1)}.$

with
$\widetilde{\mu}(s)=(\mu(s-1)^{\prime},\dots,\mu(s-p(s))^{\prime})^{\prime}$.
We can derive the closed-form solution of $\Gamma$ using that $XX^{\prime}$
can be written as

$\displaystyle
XX^{\prime}=\begin{pmatrix}\sum\limits_{n=0}^{N-1}X_{n}(1)X_{n}^{\prime}(1)&0&\dots&0\\\
0&\sum\limits_{n=0}^{N-1}X_{n}(2)X_{n}^{\prime}(2)&\ddots&\vdots\\\
\vdots&\ddots&\ddots&0\\\
0&\dots&0&\sum\limits_{n=0}^{N-1}X_{n}(S)X_{n}^{\prime}(S)\end{pmatrix}$

and that the moving average representation
$y_{Sn+s}=\mu(s)+\sum_{k=0}^{\infty}\Phi_{k}(s)\epsilon_{Sn+s-k}$ of the
PVAR($\bm{p}$) process exists. Consequently, under Assumption 3.1, we get for
$q\in\mathbb{N}$ that

$A_{q}\overset{d}{\to}\mathcal{N}(0,V_{q})\;\text{ as }\;N\to\infty,$

with

$V_{q}=\begin{pmatrix}V_{q}^{(1,1)}&V_{q}^{(1,2)}\\\
V_{q}^{(2,1)}&V^{(2,2)}\end{pmatrix}=QR_{q}\Omega_{q}R_{q}^{\prime}Q^{\prime}.$

Applying the CLT for innovations stated in Lemma A.1 to obtain $\Omega_{q}$
and using the matrices $Q$ and $R_{q}$ given above, the limiting variance
matrices $V_{q}^{(1,1)},V_{q}^{(2,1)}$ and $V^{(2,2)}$ are of dimension
$(m\sum_{s=1}^{S}(mp(s)+1)\times
m\sum_{s=1}^{S}(mp(s)+1)),(S\widetilde{m}\times m\sum_{s=1}^{S}(mp(s)+1))$ and
$(S\widetilde{m}\times S\widetilde{m})$ and given by

$\displaystyle V_{q}^{(1,1)}$
$\displaystyle=R\left[R^{\prime}\\{\Gamma\otimes\bm{I}_{m}\\}R\right]^{-1}R^{\prime}\;\bm{\Omega}_{q}^{V^{(1,1)}}\;\bigl{(}R\left[R^{\prime}\\{\Gamma\otimes\bm{I}_{m}\\}R\right]^{-1}R^{\prime}\bigr{)}^{\prime},$
$\displaystyle V_{q}^{(2,1)}$
$\displaystyle=V_{q}^{(1,2)\prime}=\bm{\Omega}_{q}^{V^{(2,1)}}\bigl{(}R\left[R^{\prime}\\{\Gamma\otimes\bm{I}_{m}\\}R\right]^{-1}R^{\prime}\bigr{)}^{\prime},$
$\displaystyle V^{(2,2)}$
$\displaystyle=L_{mS}\Bigl{(}\sum\limits_{h=-\infty}^{\infty}\tau_{0,h,0}\Bigr{)}L_{mS}^{\prime}.$

The matrix $\bm{\Omega}_{q}^{V^{(1,1)}}$ is of dimension
$(m\sum_{s=1}^{S}(mp(s)+1)\times m\sum_{s=1}^{S}(mp(s)+1))$ and can be written
as

$\bm{\Omega}_{q}^{V^{(1,1)}}=\begin{pmatrix}\bm{\Omega}_{q}^{V^{(1,1)}}(s_{1},s_{2})\\\
s_{1},s_{2}=1,\dots,S\end{pmatrix},$

where $\bm{\Omega}_{q}^{V^{(1,1)}}(s_{1},s_{2})$, $s_{1},s_{2}=1,\dots,S$ are
$(m^{2}p(s_{1})+m)\times(m^{2}p(s_{2})+m)$ dimensional and can be rewritten as

$\bm{\Omega}_{q}^{V^{(1,1)}}(s_{1},s_{2})=\begin{pmatrix}[\bm{\Omega}_{q}^{V^{(1,1)}}(s_{1},s_{2})]^{(i,j)}\\\
i,j=1,2\end{pmatrix},$

where the submatrices

$\displaystyle[\bm{\Omega}_{q}^{V^{(1,1)}}(s_{1},s_{2})]^{(1,1)}$
$\displaystyle=\Sigma_{\epsilon}(s_{1})\bm{1}(s_{1}=s_{2}),$
$\displaystyle[\bm{\Omega}_{q}^{V^{(1,1)}}(s_{1},s_{2})]^{(2,1)}$
$\displaystyle=\sum\limits_{k=1}^{q}(C_{k}(s_{1})\otimes\bm{I}_{m})\Bigl{(}\sum\limits_{h=-\infty}^{\infty}\widecheck{\kappa}_{k,h}(s_{1},s_{2})\Bigr{)},$
$\displaystyle[\bm{\Omega}_{q}^{V^{(1,1)}}(s_{1},s_{2})]^{(1,2)}$
$\displaystyle=\sum\limits_{k=1}^{q}\Bigl{(}\sum\limits_{h=-\infty}^{\infty}\widecheck{\kappa}^{\prime}_{k,h}(s_{2},s_{1})\Bigr{)}(C_{k}(s_{2})\otimes\bm{I}_{m})^{\prime},$
$\displaystyle[\bm{\Omega}_{q}^{V^{(1,1)}}(s_{1},s_{2})]^{(2,2)}$
$\displaystyle=\sum_{k,l=1}^{q}\sum_{h=-\infty}^{\infty}(C_{k}(s_{1})\otimes\bm{I}_{m})[\widecheck{\tau}_{k,h,l}(s_{1},s_{2})+\kappa_{k,h}(s_{1},s_{2})(\nu(s_{2}-l)\otimes\bm{I}_{m})^{\prime}$
$\displaystyle\quad+(\nu(s_{1}-k)\otimes\bm{I}_{m})\kappa_{l,h}(s_{1},s_{2})^{\prime}](C_{l}(s_{2})\otimes\bm{I}_{m})^{\prime}$

are of dimension $(m\times m),(m^{2}p(s_{1})\times m),(m\times m^{2}p(s_{2}))$
and $(m^{2}p(s_{1})\times m^{2}p(s_{2}))$, respectively. Notice that

$\displaystyle\widecheck{\kappa}_{k,h}(s_{1},s_{2})$
$\displaystyle=\kappa_{k,h}(s_{1},s_{2})+(\nu(s_{1}-k)\otimes\bm{I}_{m})\Sigma_{\epsilon}(s_{1})\bm{1}(h=0,s_{1}=s_{2}),$
$\displaystyle\widecheck{\tau}_{k,h,l}(s_{1},s_{2})$
$\displaystyle=\tau_{k,h,l}(s_{1},s_{2})+(\nu(s_{1}-k)\otimes\bm{I}_{m})\Sigma_{\epsilon}(s_{1})(\nu(s_{1}-l)\otimes\bm{I}_{m})^{\prime}\bm{1}(h=0,s_{1}=s_{2}).$

The $(S\widetilde{m}\times m\sum_{s=1}^{S}(mp(s)+1))$ dimensional matrix
$\bm{\Omega}_{q}^{V^{(2,1)}}$ can be decomposed into

$\bm{\Omega}_{q}^{V^{(2,1)}}=L_{mS}\bigl{(}\bigl{[}[\bm{\Omega}^{V^{(2,1)}}(s)]^{(1)}\;\;[\bm{\Omega}_{q}^{V^{(2,1)}}(s)]^{(2)}\bigr{]},\;\;s=1,\dots,S\bigr{)}$

where, for $s=1,\dots,S$, the submatrices $[\bm{\Omega}^{V^{(2,1)}}(s)]^{(1)}$
and $[\bm{\Omega}_{q}^{V^{(2,1)}}(s)]^{(2)}$ are of dimension $(Sm^{2}\times
m)$ and $(Sm^{2}\times m^{2}p(s))$ and given by

$\displaystyle[\bm{\Omega}^{V^{(2,1)}}(s)]^{(1)}$
$\displaystyle=\sum\limits_{h=-\infty}^{\infty}\begin{pmatrix}\kappa_{0,h}(1,s)\\\
\vdots\\\ \kappa_{0,h}(S,s)\end{pmatrix},$
$\displaystyle[\bm{\Omega}_{q}^{V^{(2,1)}}(s)]^{(2)}$
$\displaystyle=\sum\limits_{k=1}^{q}\sum\limits_{h=-\infty}^{\infty}\begin{pmatrix}[\tau_{0,h,k}(1,s)+\kappa_{0,h}(1,s-k)(\nu(s-k)\otimes\bm{I}_{m})^{\prime}](C_{k}(s)\otimes\bm{I}_{m})^{\prime}\\\
\vdots\\\
[\tau_{0,h,k}(S,s)+\kappa_{0,h}(S,s-k)(\nu(s-k)\otimes\bm{I}_{m})^{\prime}](C_{k}(s)\otimes\bm{I}_{m})^{\prime}\end{pmatrix}.$

Due to Assumptions 3.1 part b) follows directly since the fourth order
cumulants are assumed to be absolutely summable and the sequence of moving
average coefficient matrices $\\{\Phi_{k}(s)\\}_{k\in\mathbb{N}}$ decreases
geometrically to zero for all $s=1,\dots,S$, which implies that
$\\{C_{k}(s)\\}_{k\in\mathbb{N}}$, $s=1,\dots,S$ decreases geometrically to
zero as well. Hence, we have

$V_{q}\to V=\begin{pmatrix}V^{(1,1)}&V^{(1,2)}\\\
V^{(2,1)}&V^{(2,2)}\end{pmatrix}\;\text{ as }\;q\to\infty.$

Note that the sums $\sum_{k=1}^{q}$ and $\sum_{k,l=1}^{q}$ in $V_{q}^{(1,1)}$,
$V_{q}^{(1,2)}$ and $V_{q}^{(2,1)}$ are replaced by $\sum_{k=1}^{\infty}$ and
$\sum_{k,l=1}^{\infty}$, respectively, to obtain the limiting variances
$V^{(1,1)}$, $V^{(2,1)}$ and $V^{(1,2)}$. The variance matrix $V^{(2,2)}$ is
provided above. To prove c), note that odd parts of the first subvector and
the complete second subvector of $A-A_{q}$ in (A.2) are zero. Therefore, it
suffices to show part c) only for the even parts of the first subvector of
$A-A_{q}$. For any constant vector $d\in\mathbb{R}^{m^{2}p(s)}$ and
$\delta>0$, we use the Markov inequality to get

$\displaystyle
P\Biggl{(}\Biggl{|}\sum\limits_{k=q+1}^{\infty}d^{\prime}(C_{k}(s)\otimes\bm{I}_{m})\frac{1}{\sqrt{N}}\sum\limits_{n=0}^{N-1}vec\\{\epsilon_{Sn+s}[\epsilon_{Sn+s-k}^{\prime}+\nu(s-k)^{\prime}]\\}\Biggr{|}>\delta\Biggr{)}$
$\displaystyle\leq\;\;$
$\displaystyle\frac{1}{\delta^{2}N}E\Biggl{(}\Biggl{|}\sum\limits_{k=q+1}^{\infty}d^{\prime}(C_{k}(s)\otimes\bm{I}_{m})\sum\limits_{n=0}^{N-1}vec\\{\epsilon_{Sn+s}[\epsilon_{Sn+s-k}^{\prime}+\nu(s-k)^{\prime}]\\}\Biggr{|}^{2}\Biggr{)}$
$\displaystyle=\;\;$
$\displaystyle\frac{1}{\delta^{2}}\sum\limits_{k,l=q+1}^{\infty}d^{\prime}(C_{k}(s)\otimes\bm{I}_{m})\Biggl{(}\frac{1}{N}\sum\limits_{n_{1},n_{2}=0}^{N-1}E\Bigl{(}vec\\{\epsilon_{Sn+s}[\epsilon_{Sn+s-k}^{\prime}+\nu(s-k)^{\prime}]\\}$
$\displaystyle\times
vec\\{\epsilon_{Sn+s}[\epsilon_{Sn+s-k}^{\prime}+\nu(s-k)^{\prime}]\\}^{\prime}\Bigr{)}\Biggr{)}(C_{l}(s)\otimes\bm{I}_{m})^{\prime}d$
$\displaystyle\underset{N\to\infty}{\to}$
$\displaystyle\;\;\frac{1}{\delta^{2}}\sum\limits_{k,l=q+1}^{\infty}d^{\prime}(C_{k}(s)\otimes\bm{I}_{m})\Biggl{(}\sum\limits_{h=-\infty}^{\infty}\Bigl{[}\widecheck{\tau}_{k,h,l}(s,s)+\kappa_{k,h}(s,s)(\nu(s-l)\otimes\bm{I}_{m})^{\prime}$
$\displaystyle+(\nu(s-k)\otimes\bm{I}_{m})\kappa_{l,h}(s,s)^{\prime}\Bigr{]}\Biggr{)}(C_{l}(s)\otimes\bm{I}_{m})^{\prime}d$
$\displaystyle\underset{q\to\infty}{\to}$ $\displaystyle\;\;0$

for all $s=1,\dots,S$. Hence, we proved

$\sqrt{N}\begin{pmatrix}\widehat{\beta}_{res}-\beta\\\
\widetilde{\sigma}-\sigma\end{pmatrix}\overset{d}{\to}\mathcal{N}\Biggl{(}0,\begin{pmatrix}V^{(1,1)}&V^{(1,2)}\\\
V^{(2,1)}&V^{(2,2)}\end{pmatrix}\Biggr{)}$

which, for fixed $S$, implies

$\sqrt{T}\begin{pmatrix}\widehat{\beta}_{res}-\beta\\\
\widetilde{\sigma}-\sigma\end{pmatrix}\overset{d}{\to}\mathcal{N}\Biggl{(}0,S\begin{pmatrix}V^{(1,1)}&V^{(1,2)}\\\
V^{(2,1)}&V^{(2,2)}\end{pmatrix}\Biggr{)}.$

######  Lemma A.1 (CLT for Innovations).

Let $q\in\mathbb{N}$ and define
$W_{n,q}=(W^{(1)\prime}_{n,q},W^{(2)\prime}_{n})^{\prime}$ as

$\displaystyle W^{(1)}_{n,q}=\begin{pmatrix}W^{(1,1)}_{n,q}\\\ \vdots\\\
W^{(1,S)}_{n,q}\end{pmatrix},\quad
W^{(2)\prime}_{n}=L_{mS}\begin{pmatrix}vec\\{\epsilon_{Sn+1}\epsilon_{Sn+1}^{\prime}\\}-vec\\{\Sigma_{\epsilon}(1)\\}\\\
\vdots\\\
vec\\{\epsilon_{Sn+S}\epsilon_{Sn+S}^{\prime}\\}-vec\\{\Sigma_{\epsilon}(S)\\}\end{pmatrix},$

where, for $s=1,\dots,S$,

$\displaystyle W^{(1,s)}_{n,q}=\begin{pmatrix}\epsilon_{Sn+s}\\\
vec\\{\epsilon_{Sn+s}[\epsilon_{Sn+s-1}^{\prime}+\nu(s-1)^{\prime}]\\}\\\
\vdots\\\
vec\\{\epsilon_{Sn+s}[\epsilon_{Sn+s-q}^{\prime}+\nu(s-q)^{\prime}]\\}\end{pmatrix}.$

Then, under Assumption 3.1, we get the following CLT

$\frac{1}{\sqrt{N}}\sum_{n=0}^{N-1}\begin{pmatrix}W_{n,q}^{(1)}\\\
W_{n,q}^{(2)}\end{pmatrix}\overset{d}{\to}\mathcal{N}\bigl{(}\bm{0},\bm{\Omega}_{q}\bigr{)},$

with
$(S(m+qm^{2})+S\widetilde{m})\times(S(m+qm^{2})+S\widetilde{m})$-dimensional
limiting variance matrix $\bm{\Omega}_{q}$.

###### Proof.

To prove Lemma A.1, we state the limiting variance matrix $\bm{\Omega}_{q}$ as

$\bm{\Omega}_{q}=\begin{pmatrix}\bm{\Omega}_{q}^{(1)}&\bm{\Omega}_{q}^{(2,1)\prime}\\\
\bm{\Omega}_{q}^{(2,1)}&\bm{\Omega}^{(2)}\end{pmatrix},$

where the matrices $\bm{\Omega}_{q}^{(1)},\bm{\Omega}_{q}^{(2,1)}$ and
$\bm{\Omega}^{(2)}$ are of dimension $(S(m+qm^{2})\times
S(m+qm^{2})),(S\widetilde{m}\times S(m+qm^{2}))$ and $(S\widetilde{m}\times
S\widetilde{m})$. The matrix $\bm{\Omega}_{q}^{(1)}$ can be represented as

$\displaystyle\bm{\Omega}_{q}^{(1)}=\begin{pmatrix}\bm{\Omega}_{q}^{(1)}(s_{1},s_{2})\\\
s_{1},s_{2}=1,\dots,S\end{pmatrix},$

where $\bm{\Omega}_{q}^{(1)}(s_{1},s_{2})$ can be decomposed as

$\displaystyle\bm{\Omega}_{q}^{(1)}(s_{1},s_{2})=\begin{pmatrix}\Xi_{1}^{(1,1)}(s_{1},s_{2})&\Xi_{1,q}^{(1,2)}(s_{1},s_{2})\\\
\Xi_{1,q}^{(2,1)}(s_{1},s_{2})&\Xi_{1,q}^{(2,2)}(s_{1},s_{2})\end{pmatrix}\in\mathbb{R}^{(m+qm^{2})\times(m+qm^{2})},\;\;\;s_{1},s_{2}=1,\dots,S.$

For $s_{1},s_{2}=1,\dots,S$, the submatrices
$\Xi_{1}^{(1,1)}(s_{1},s_{2}),\Xi_{1,q}^{(2,1)}(s_{1},s_{2}),\Xi_{1,q}^{(1,2)}(s_{1},s_{2})$
and $\Xi_{1,q}^{(2,2)}(s_{1},s_{2})$ are of dimension $(m\times
m),(qm^{2}\times m)$, $(m\times qm^{2})$ and $(qm^{2}\times qm^{2})$ and given
by

$\displaystyle\Xi_{1}^{(1,1)}(s_{1},s_{2})$
$\displaystyle=\sum_{h=-\infty}^{\infty}Cov(\epsilon_{Sn+s_{1}},\epsilon_{S(n-h)+s_{2}})=\Sigma_{\epsilon}(s_{1})\bm{1}(s_{1}=s_{2}),$
$\displaystyle\Xi_{1,q}^{(2,1)}(s_{1},s_{2})$
$\displaystyle=\sum\limits_{h=-\infty}^{\infty}\begin{pmatrix}Cov\bigl{(}vec\\{\epsilon_{Sn+s_{1}}[\epsilon_{Sn+s_{1}-1}^{\prime}+\nu(s_{1}-1)^{\prime}],\epsilon_{S(n-h)+s_{2}}\\}\bigr{)}\\\
\vdots\\\
Cov\bigl{(}vec\\{\epsilon_{Sn+s_{1}}[\epsilon_{Sn+s_{1}-q}^{\prime}+\nu(s_{1}-q)^{\prime}]\\},\epsilon_{S(n-h)+s_{2}}\bigr{)}\end{pmatrix}$
$\displaystyle=\sum\limits_{h=-\infty}^{\infty}\begin{pmatrix}\kappa_{1,h}(s_{1},s_{2})+(\nu(s_{1}-1)\otimes\bm{I}_{m})\Sigma_{\epsilon}(s_{1})\bm{1}(h=0,s_{1}=s_{2})\\\
\vdots\\\
\kappa_{q,h}(s_{1},s_{2})+(\nu(s_{1}-q)\otimes\bm{I}_{m})\Sigma_{\epsilon}(s_{1})\bm{1}(h=0,s_{1}=s_{2})\end{pmatrix},$
$\displaystyle\Xi_{1,q}^{(1,2)}(s_{1},s_{2})$
$\displaystyle=\sum\limits_{h=-\infty}^{\infty}\begin{pmatrix}Cov\bigl{(}vec\\{\epsilon_{Sn+s_{2}}[\epsilon_{Sn+s_{2}-1}^{\prime}+\nu(s_{2}-1)^{\prime}],\epsilon_{S(n-h)+s_{1}}\\}\bigr{)}\\\
\vdots\\\
Cov\bigl{(}vec\\{\epsilon_{Sn+s_{2}}[\epsilon_{Sn+s_{2}-q}^{\prime}+\nu(s_{2}-q)^{\prime}]\\},\epsilon_{S(n-h)+s_{1}}\bigr{)}\end{pmatrix}^{\prime}$
$\displaystyle=\sum\limits_{h=-\infty}^{\infty}\begin{pmatrix}\kappa_{1,h}(s_{2},s_{1})+(\nu(s_{1}-1)\otimes\bm{I}_{m})\Sigma_{\epsilon}(s_{1})\bm{1}(h=0,s_{1}=s_{2})\\\
\vdots\\\
\kappa_{q,h}(s_{2},s_{1})+(\nu(s_{1}-q)\otimes\bm{I}_{m})\Sigma_{\epsilon}(s_{1})\bm{1}(h=0,s_{1}=s_{2})\end{pmatrix}^{\prime},$
$\displaystyle\Xi_{1,q}^{(2,2)}(s_{1},s_{2})$
$\displaystyle=\sum\limits_{h=-\infty}^{\infty}\begin{pmatrix}Cov\bigl{(}vec\\{\epsilon_{Sn+s_{1}}[\epsilon_{Sn+s_{1}-k}^{\prime}+\nu(s_{1}-k)^{\prime}]\\},\\\
vec\\{\epsilon_{S(n-h)+s_{2}}[\epsilon_{S(n-h)+s_{2}-l}^{\prime}+\nu(s_{2}-l)^{\prime}]\\}\bigr{)}\\\
k,l=1,\dots,q\end{pmatrix}$
$\displaystyle=\sum\limits_{h=-\infty}^{\infty}\begin{pmatrix}\tau_{k,h,l}(s_{1},s_{2})+\kappa_{k,h}(s_{1},s_{2})(\nu(s_{2}-l)\otimes\bm{I}_{m})^{\prime}+(\nu(s_{1}-k)\otimes\bm{I}_{m})\kappa_{l,h}(s_{1},s_{2})^{\prime}\\\
+(\nu(s_{1}-k)\otimes\bm{I}_{m})\Sigma_{\epsilon}(s_{1})(\nu(s_{1}-l)\otimes\bm{I}_{m})^{\prime}\bm{1}(h=0,s_{1}=s_{2})\\\
k,l=1,\dots,q\end{pmatrix}.$

The matrix $\bm{\Omega}^{(2)}$ is given by

$\displaystyle\bm{\Omega}^{(2)}$
$\displaystyle=L_{mS}\left(\sum_{h=-\infty}^{\infty}Cov\left(\begin{pmatrix}vec\\{\epsilon_{Sn+1}\epsilon_{Sn+1}^{\prime}\\}\\\
\vdots\\\
vec\\{\epsilon_{Sn+S}\epsilon_{Sn+S}^{\prime}\\}\end{pmatrix},\begin{pmatrix}vec\\{\epsilon_{S(n-h)+1}\epsilon_{S(n-h)+1}^{\prime}\\}\\\
\vdots\\\
vec\\{\epsilon_{S(n-h)+S}\epsilon_{S(n-h)+S}^{\prime}\\}\end{pmatrix}\right)\right)L_{mS}^{\prime}$
$\displaystyle=L_{mS}\Bigl{(}\sum\limits_{h=-\infty}^{\infty}\tau_{0,h,0}\Bigr{)}L_{mS}^{\prime},\;\;\;\in\mathbb{R}^{S\tilde{m}\times
S\tilde{m}}.$

The $(S\widetilde{m}\times S(m+qm^{2}))$-dimensional matrix
$\bm{\Omega}_{q}^{(2,1)}$ can be represented by

$\bm{\Omega}_{q}^{(2,1)}=L_{mS}\left(\left[\Xi_{2,1}^{(1)}(s)\;\;\Xi_{2,1,q}^{(2)}(s)\right],\;\;s=1,\dots,S\right),$

where, for $s=1,\dots,S$, the submatrices $\Xi_{2,1}^{(1)}(s)$ and
$\Xi_{2,1,q}^{(2)}(s)$ are of dimension $(Sm^{2}\times m)$ and $(Sm^{2}\times
qm^{2})$ and given by

$\displaystyle\Xi_{2,1}^{(1)}(s)$
$\displaystyle=\sum\limits_{h=-\infty}^{\infty}Cov\Biggl{(}\begin{pmatrix}vec\\{\epsilon_{Sn+1}\epsilon_{Sn+1}^{\prime}\\}\\\
\vdots\\\
vec\\{\epsilon_{Sn+S}\epsilon_{Sn+S}^{\prime}\\}\end{pmatrix},\epsilon_{S(n-h)+s}\Biggr{)}=\sum\limits_{h=-\infty}^{\infty}\begin{pmatrix}\kappa_{0,h}(1,s)\\\
\vdots\\\ \kappa_{0,h}(S,s)\end{pmatrix}$ $\displaystyle\Xi_{2,1,q}^{(2)}(s)$
$\displaystyle=\sum\limits_{h=-\infty}^{\infty}Cov\Biggl{(}\begin{pmatrix}vec\\{\epsilon_{Sn+1}\epsilon_{Sn+1}^{\prime}\\}\\\
\vdots\\\
vec\\{\epsilon_{Sn+S}\epsilon_{Sn+S}^{\prime}\\}\end{pmatrix},\begin{pmatrix}vec\\{\epsilon_{S(n-h)+s}[\epsilon_{S(n-h)+s-1}^{\prime}+\nu(s-1)^{\prime}]\\}\\\
\vdots\\\
vec\\{\epsilon_{S(n-h)+s}[\epsilon_{S(n-h)+s-q}^{\prime}+\nu(s-q)^{\prime}]\\}\end{pmatrix}\Biggr{)}$
$\displaystyle=\sum\limits_{h=-\infty}^{\infty}\begin{pmatrix}\tau_{0,h,k}(1,s)+\kappa_{0,h}(1,s-k)(\nu(s-k)\otimes\bm{I}_{m})^{\prime}&\\\
\vdots&,\;\;\;k=1,\dots,q\\\
\tau_{0,h,k}(S,s)+\kappa_{0,h}(S,s-k)(\nu(s-k)\otimes\bm{I}_{m})^{\prime}&\end{pmatrix}.$

∎

##  Appendix B Joint CLT with Periodic Strong White Noise

######  Corollary B.1 (Joint CLT with Periodic Strong White Noise).

Under Assumption 3.1 with $(ii)$ and $(iii)$ replaced by $(ii)^{\prime}$ and
$(iii)^{\prime}$, we have

$\sqrt{T}\begin{pmatrix}\widehat{\beta}_{res}-\beta\\\
\widehat{\sigma}-\sigma\end{pmatrix}\overset{d}{\to}\mathcal{N}(\bm{0},SV_{iid}),\quad
V_{iid}=\begin{pmatrix}V^{(1,1)}_{iid}&V^{(1,2)}_{iid}\\\
V^{(2,1)}_{iid}&V^{(2,2)}_{iid}\end{pmatrix},$

where, using $\frac{1}{N}XX^{\prime}\overset{p}{\to}\Gamma$, the submatrices
$V^{(1,1)}_{iid},V^{(2,1)}_{iid},V^{(1,2)}_{iid}$ and $V^{(2,2)}_{iid}$ are
given by

$\displaystyle V^{(1,1)}_{iid}$
$\displaystyle=R\left[R^{\prime}\\{\Gamma\otimes\bm{I}_{m}\\}R\right]^{-1}R^{\prime}\;\bm{\Omega}^{V^{(1,1)}}_{iid}\bigl{(}R\left[R^{\prime}\\{\Gamma\otimes\bm{I}_{m}\\}R\right]^{-1}R^{\prime}\bigr{)}^{\prime},$
$\displaystyle V^{(2,1)}_{iid}$
$\displaystyle=V^{(1,2)\prime}_{iid}=\bm{0}_{S\widetilde{m}\times
m\sum_{s=1}^{S}(mp(s)+1)},$ $\displaystyle V^{(2,2)}_{iid}$
$\displaystyle=L_{mS}\\{\tau_{0,0,0}\\}L_{mS}^{\prime},$

where

$\bm{\Omega}^{V^{(1,1)}}_{iid}=\begin{pmatrix}[\bm{\Omega}^{V^{(1,1)}}_{iid}(1)]^{(1,1)}&[\bm{\Omega}^{V^{(1,1)}}_{iid}(1)]^{(2,1)\prime}&\bm{0}&\dots&\bm{0}\\\
[\bm{\Omega}^{V^{(1,1)}}_{iid}(1)]^{(2,1)}&[\bm{\Omega}^{V^{(1,1)}}_{iid}(1)]^{(2,2)}&&\ddots&\vdots\\\
\bm{0}&&\ddots&&\bm{0}\\\
\vdots&\ddots&&[\bm{\Omega}^{V^{(1,1)}}_{iid}(S)]^{(1,1)}&[\bm{\Omega}^{V^{(1,1)}}_{iid}(S)]^{(2,1)\prime}\\\
\bm{0}&\dots&\bm{0}&[\bm{\Omega}^{V^{(1,1)}}_{iid}(S)]^{(2,1)}&[\bm{\Omega}^{V^{(1,1)}}_{iid}(S)]^{(2,2)}\end{pmatrix}$

with, for $s=1,\dots,S$,

$\displaystyle[\bm{\Omega}^{V^{(1,1)}}_{iid}(s)]^{(1,1)}$
$\displaystyle=\Sigma_{\epsilon}(s),$
$\displaystyle[\bm{\Omega}^{V^{(1,1)}}_{iid}(s)]^{(2,1)}$
$\displaystyle=\sum\limits_{i=1}^{\infty}(C_{i}(s)\otimes\bm{I}_{m})[\nu(s-i)\otimes\Sigma_{\epsilon}(s)],$
$\displaystyle[\bm{\Omega}^{V^{(1,1)}}_{iid}(s)]^{(2,2)}$
$\displaystyle=\sum\limits_{i=1}^{\infty}(C_{i}(s)\otimes\bm{I}_{m})[\tau_{i,0,i}(s,s)+[\nu(s-i)\nu(s-i)^{\prime}\otimes\Sigma_{\epsilon}(s)]](C_{i}(s)\otimes\bm{I}_{m})^{\prime}.$

In contrast to the case with periodic weak white noise addressed in Theorem
3.1, in Corollary B.1, the off-diagonal sub-matrices
$V^{(2,1)}_{iid}=V^{(1,2)\prime}_{iid}$ are zero-matrices leading to
asymptotically independent estimators $\widehat{\beta}_{res}$ and
$\widehat{\sigma}$. Since $\Gamma$ is block diagonal, we obtain the same
results as Ursu & Duchesne (2009), including asymptotic independence among the
seasonal PVAR estimators $\widehat{\beta}_{res}(s)$, $s=1,\dots,S$, if no
restrictions _across_ seasons are imposed. Further, the matrix $\tau_{0,0,0}$
is also block diagonal, since $\tau_{0,0,0}(s_{1},s_{2})=\bm{0}_{m^{2}\times
m^{2}}$ for $s_{1}\neq s_{2}$ implying also asymptotic independence among the
seasonal covariance estimators $\widehat{\sigma}(s),s=1,\dots,S$.

##  Appendix C Proof of Theorem 4.1

Theorem 4.1 states that under Assumptions 3.1, we have

$\displaystyle\sqrt{T}\bigl{(}vec\\{\widehat{\Pi}_{h}^{IR}-\Pi_{h}^{IR}\\}\bigr{)}\overset{d}{\to}\mathcal{N}\bigl{(}\bm{0},S\;\Sigma_{h}^{\Pi}\bigr{)},$
$\displaystyle\sqrt{T}\bigl{(}vec\\{\widehat{\Psi}_{h}^{SIR}-\Psi_{h}^{SIR}\\}\bigr{)}\overset{d}{\to}\mathcal{N}\bigl{(}\bm{0},S\;\Sigma_{h}^{\Psi}\bigr{)},\;\;\;h\in\mathbb{N}_{0},$

where

$\displaystyle\Sigma_{h}^{\Pi}$
$\displaystyle=G^{\Pi}_{h}V^{(1,1)}G_{h}^{\Pi\prime},$
$\displaystyle\Sigma_{h}^{\Psi}$
$\displaystyle=F_{h}V^{(1,1)}F_{h}^{\prime}+D_{h}V^{(2,2)}D_{h}^{\prime}+F_{h}V^{(2,1)\prime}D_{h}^{\prime}+D_{h}V^{(2,1)}F_{h}^{\prime}$

with $G_{h}=\frac{\partial vec\\{\Pi_{h}^{IR}\\}}{\partial\beta^{\prime}}$,
$F_{h}=\frac{\partial vec\\{\Psi_{h}^{SIR}\\}}{\partial\beta^{\prime}}$ and
$D_{h}=\frac{\partial vec\\{\Psi_{h}^{SIR}\\}}{\partial\sigma^{\prime}}$.

###### Proof.

In order to prove the CLT of the (structural) impulse responses stated in
Theorem 4.1, we make use of Proposition 3.6 of Lütkepohl (2005) and Theorem
3.1. For $h\in\mathbb{N}_{0}$, from the relation
$\Pi_{h}^{IR}=\Pi_{h}\bm{A}_{0}^{-1}$, we immediately get that the first
derivative of $\Pi_{h}^{IR}$ with respect to $\beta$ is given by

$\displaystyle G_{h}=\frac{\partial
vec\\{\Pi_{h}^{IR}\\}}{\partial\beta^{\prime}}=\frac{\partial
vec\\{\Pi_{h}\bm{A}_{0}^{-1}\\}}{\partial\beta^{\prime}}=(\bm{I}_{Sm}\otimes\Pi_{h})\frac{\partial
vec\\{\bm{A}_{0}^{-1}\\}}{\partial\beta^{\prime}}+(\bm{A}_{0}^{-1\prime}\otimes\bm{I}_{Sm})\frac{\partial
vec\\{\Pi_{h}\\}}{\partial\beta^{\prime}},$

where

$\displaystyle\frac{\partial
vec\\{\bm{A}_{0}^{-1}\\}}{\partial\beta^{\prime}}=-(\bm{A}_{0}^{-1\prime}\otimes\bm{A}_{0}^{-1})\frac{\partial
vec\\{\bm{A}_{0}\\}}{\partial\beta^{\prime}}=-(\bm{A}_{0}^{-1\prime}\otimes\bm{A}_{0}^{-1})J_{\bm{A}_{0}}$

and

$\displaystyle\frac{\partial vec\\{\Pi_{h}\\}}{\partial\beta^{\prime}}$
$\displaystyle=\frac{\partial
vec\\{J\bm{A}^{h}J^{\prime}\\}}{\partial\beta^{\prime}}=(J\otimes
J)\frac{\partial vec\\{\bm{A}^{h}\\}}{\partial\beta^{\prime}}=(J\otimes
J)\Bigl{[}\sum\limits_{i=0}^{h-1}(\bm{A}^{\prime})^{h-1-i}\otimes\bm{A}^{h}\Bigr{]}\frac{\partial
vec\\{\bm{A}\\}}{\partial\beta^{\prime}}$
$\displaystyle=\Bigl{[}\sum\limits_{i=0}^{h-1}J(\bm{A}^{\prime})^{h-1-i}\otimes
J\bm{A}^{h}\Bigr{]}J_{\bm{A}}.$

Here, we define $J_{\bm{A}_{0}}=\frac{\partial
vec\\{\bm{A}_{0}\\}}{\partial\beta^{\prime}}$ and $J_{\bm{A}}=\frac{\partial
vec\\{\bm{A}\\}}{\partial\beta^{\prime}}$. Note that $\Pi_{h}$ can be
represented as $\Pi_{h}=J\bm{A}^{h}J^{\prime}$ with
$J=[\bm{I}_{Sm},\bm{0}_{Sm},\dots,\bm{0}_{Sm}]$ and VAR($P$) companion matrix
$\bm{A}$. We can represent the companion matrix $\bm{A}$ of the time-invariant
VAR($P$) as the product $\bm{A}=\bm{A}_{0}^{*}\bm{A}^{*}$ where

$\displaystyle\bm{A}_{0}^{*}=\begin{pmatrix}\bm{A}_{0}^{-1}&\bm{0}&\dots&\bm{0}\\\
\bm{0}&\bm{I}_{Sm}&\ddots&\vdots\\\ \vdots&\ddots&\ddots&\bm{0}\\\
\bm{0}&\dots&\bm{0}&\bm{I}_{Sm}\end{pmatrix}\quad\text{and}\quad\bm{A}^{*}=\begin{pmatrix}\bm{A}_{1}&\bm{A}_{2}&\dots&\bm{A}_{P}\\\
\bm{I}_{Sm}&\bm{0}&\dots&\bm{0}\\\ &\ddots&\ddots&\vdots\\\
\bm{0}&&\bm{I}_{Sm}&\bm{0}\end{pmatrix}$

and use this representation in order to calculate $J_{\bm{A}}$. Hence, $G_{h}$
can be represented as

$G_{h}=-(\bm{I}_{Sm}\otimes\Pi_{h})(\bm{A}_{0}^{-1\prime}\otimes\bm{A}_{0}^{-1})J_{\bm{A}_{0}}+(\bm{A}_{0}^{-1\prime}\otimes\bm{I}_{Sm})\Bigl{[}\sum\limits_{i=0}^{h-1}J(\bm{A}^{\prime})^{h-1-i}\otimes
J\bm{A}^{h}\Bigr{]}J_{\bm{A}}.$

From $\Psi_{h}^{SIR}=\Pi_{h}^{IR}H_{0}$, we immediately get

$\displaystyle F_{h}$ $\displaystyle=\frac{\partial
vec\\{\Psi_{h}^{SIR}\\}}{\partial\beta^{\prime}}=\frac{\partial
vec\\{\Pi_{h}^{IR}H_{0}\\}}{\partial\beta^{\prime}}=(\bm{I}_{Sm}\otimes\Pi_{h}^{IR})\frac{\partial
vec\\{H_{0}\\}}{\partial\beta^{\prime}}+(H_{0}^{\prime}\otimes\bm{I}_{Sm})G_{h}$
$\displaystyle=(H_{0}^{\prime}\otimes\bm{I}_{Sm})G_{h}.$