diff --git a/.gitignore b/.gitignore
index eee0010..e81e5c5 100644
--- a/.gitignore
+++ b/.gitignore
@@ -316,3 +316,5 @@ TSWLatexianTemp*
 *.glstex
 
 # End of https://www.gitignore.io/api/latex
+tests/test_file_iv.R
+tests/iv_testfile.R
diff --git a/DESCRIPTION b/DESCRIPTION
index 5e4c333..93d38a9 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -6,12 +6,15 @@ Date: 2024-02-13
 Authors@R: c(person("Nikolas", "Kuschnig", role = c("aut", "cre"), email = "nikolas.kuschnig@wu.ac.at", comment = c(ORCID = "0000-0002-6642-2543")),
     person("Lukas", "Vashold", role = "aut", comment = c(ORCID = "0000-0002-3562-3414")),
     person("Nirai", "Tomass", role = "ctb"),
+    person("Gabriel", "Konecny", role = "ctb"),
     person("Michael", "McCracken", role = "dtc"),
     person("Serena", "Ng", role = "dtc"))
 Author: Nikolas Kuschnig [aut, cre] (<https://orcid.org/0000-0002-6642-2543>),
     Lukas Vashold [aut] (<https://orcid.org/0000-0002-3562-3414>),
     Nirai Tomass [ctb],
-    Michael McCracken [dtc], Serena Ng [dtc]
+    Gabriel Konecny [ctb],
+    Michael McCracken [dtc], Serena Ng [dtc],
+    Silvia Miranda-Agrippino [dtc], Giovanni Ricco [dtc], Riccardo Degasperi [dtc] 
 Maintainer: Nikolas Kuschnig <nikolas.kuschnig@wu.ac.at>
 Description: Estimation of hierarchical Bayesian vector autoregressive models
     following Kuschnig & Vashold (2021) <doi:10.18637/jss.v100.i14>.
@@ -28,4 +31,4 @@ Suggests: coda, vars, tinytest
 License: GPL-3 | file LICENSE
 Encoding: UTF-8
 LazyData: true
-RoxygenNote: 7.3.1
+RoxygenNote: 7.3.2
diff --git a/NAMESPACE b/NAMESPACE
index b11787b..14cb865 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -50,6 +50,7 @@ S3method(print,bvar_irf_summary)
 S3method(print,bvar_resid)
 S3method(print,bvar_summary)
 S3method(print,bvar_vcovs)
+S3method(print,check_iv)
 S3method(residuals,bvar)
 S3method(rmse,bvar)
 S3method(rmse,default)
diff --git a/R/10_bvar.R b/R/10_bvar.R
index 3149f68..fedb578 100644
--- a/R/10_bvar.R
+++ b/R/10_bvar.R
@@ -187,6 +187,10 @@ bvar <- function(
   variables <- name_deps(variables = colnames(data), M = M)
   explanatories <- name_expl(variables = variables, M = M, lags = lags)
 
+  # if(!is.null(irf$instrument)){instrument <- irf$instrument
+  # check_iv_result <- check_iv(Y, instrument)
+  # print(check_iv_result)}
+
 
   # Priors -----
 
diff --git a/R/12_aux.R b/R/12_aux.R
index 75469c3..38a12ff 100644
--- a/R/12_aux.R
+++ b/R/12_aux.R
@@ -267,3 +267,124 @@ quantile_check <- function(conf_bands) {
 
   return(quants)
 }
+
+#' Set dates
+#'
+#' Changes rownames to a character string containing dates. Common format of
+#' dates or other index is needed if the length of the \emph{instrument} and
+#' reduced form residuals differ. In such case, the identification is based on
+#' a common subset of residuals and the instrument.
+#'
+#' @param data Numeric vector, matrix or a data.frame. Note that observations
+#' are expected to be ordered from earliest to latest, and variables in the
+#' columns.
+#' @param start Starting date of the data in form "YYYY-MM-DD".
+#' @param frequency Frequency of the data: "year", "quarter", "month", "week" or "day".
+#'
+#' @return Returns data with rownames containing sequence of dates as character.
+set_dates <- function(data, start, frequency) {
+
+  # Check if the input is a data frame, matrix, or numeric vector
+  if (!is.data.frame(data) && !is.matrix(data) && !(is.vector(data) && is.numeric(data))) {
+    stop("Error: data must be either a data frame, a matrix, or a numeric vector.")
+  }
+
+  # Determine the number of elements to set row names for
+  num_elements <- if (is.data.frame(data)) {
+    nrow(data)
+  } else if (is.matrix(data)) {
+    nrow(data)
+  } else {
+    length(data)
+  }
+
+  # Set the row names or names based on the sequence of dates
+  rownames(data) <- seq(as.Date(start), by = frequency, length.out = num_elements)
+
+  return(data)
+}
+
+
+#' Check IV
+#'
+#' Intersects the \emph{residuals} and the \emph{instrument} based on an index
+#' (e.g. dates) provided in rownames and names, respectively. Returns an error
+#'  if the intersection is empty.
+#'
+#' @param residuals Data frame or numeric matrix. During estimation, draw specific residuals
+#' are used and rownames are inherited from data.
+#' @param instrument Numeric vector.
+#' @param proxyvar character string. Variable for which instrument is provided -
+#' Needs to correspond to one of the column names of the data used in
+#' \emph{bvar}.
+#' @param manual_matching If set to TRUE, user is not expected to specify common
+#' index for the data and instrument. Instead the length of instrument needs to
+#' match the length of residuals (= length(data inputted in bvar) - lags).
+#' @return Returns a list of class \code{check_iv} with information on
+#' values and indices corresponding to the intersection. Additionally the list
+#' contains inputs provided to the function.
+#' @noRd
+check_iv <- function(residuals, instrument, manual_matching = FALSE,
+                     proxyvar) {
+
+residuals <- data.frame(residuals)
+
+if(manual_matching == TRUE){
+  dates_residuals <- NULL
+  dates_instrument <- NULL
+  common_dates <- NULL
+
+  if(!nrow(residuals) == length(instrument)){
+  stop("manual_matching is TRUE but the length of residuals does not match the length of instrument!
+  Note that residuals are shorter than the data inputted to bvar, since first l observations of data are not part of residuals,
+  where l corresponds to number of lags used.")
+  }
+
+  residuals_shortened <- residuals
+  instrument_shortened <- instrument
+} else{
+
+  if (!proxyvar %in% colnames(residuals)) {
+    stop(paste("The proxied variable", proxyvar, "was not found in the column names of residuals."))
+  }
+  # Get the dates or other indices
+  dates_residuals <- rownames(residuals)
+  dates_instrument <- names(instrument)
+
+  # Find the intersection of dates (common dates)
+  common_dates <- intersect(dates_residuals, dates_instrument)
+
+  # If there are no matches and throw an error
+  if (length(common_dates) == 0) {
+    stop("Could not match instrument with the data; please provide common indices for both. See `bv_irf` for details.")
+    }
+
+  # Shorten both objects to the common dates
+  residuals_shortened <- residuals[common_dates,]
+  instrument_shortened <- instrument[common_dates]
+}
+
+  out <- structure(mget(c("residuals", "instrument",
+                          "dates_residuals", "dates_instrument",
+                          "residuals_shortened", "instrument_shortened",
+                          "common_dates", "manual_matching")),
+                       class = "check_iv")
+
+  return(out)
+}
+
+
+get_swapped_index <- function(residuals, proxyvar) {
+  # Find the column index of the proxy variable
+  j <- which(colnames(residuals) == proxyvar)
+
+  # Create a sequence of column indices
+  index <- seq_len(ncol(residuals))
+
+  # Swap the first column index with the proxyvar column index
+  index[1] <- index[j]
+  index[j] <- 1
+
+  return(index)
+}
+
diff --git a/R/60_irf_setup.R b/R/60_irf_setup.R
index f91aae1..8863633 100644
--- a/R/60_irf_setup.R
+++ b/R/60_irf_setup.R
@@ -8,15 +8,25 @@
 #' as well as if and what kind of identification should be used. See the Details
 #' section for further information on identification. Identification can be
 #' achieved via Cholesky decomposition, sign restrictions (Rubio-Ramirez,
-#' Waggoner and Zha, 2010), and zero and sign restrictions (Arias,
-#' Rubio-Ramirez and Waggoner, 2018).
+#' Waggoner and Zha, 2010), combination of zero and sign restrictions (Arias,
+#' Rubio-Ramirez and Waggoner, 2018) or using an external instrument
+#' (Stock and Watson 2012, Mertens and Ravn 2013)
 #'
 #' Identification can be performed via Cholesky decomposition, sign
-#' restrictions, or zero and sign restrictions. The algorithm
-#' for generating suitable sign restrictions follows Rubio-Ramirez, Waggoner
-#' and Zha (2010), while the one for zero and sign restrictions follows
-#' Arias, Rubio-Ramirez and Waggoner (2018).
+#' restrictions, zero and sign restrictions or using an external instrument.
+#' The algorithm for generating suitable sign restrictions follows
+#' Rubio-Ramirez, Waggoner and Zha (2010), while the one for zero and sign
+#' restrictions follows Arias, Rubio-Ramirez and Waggoner (2018).
 #' Note the possiblity of finding no suitable zero/sign restrictions.
+#' The proxy SVAR identification follows the implementation in
+#' Miranda-Agrippino and Ricco (2021) and examples make use of their instrument
+#' extended in Degasperi and Ricco (2021). The implemented proxy SVAR allows
+#' for identification using subsample of residuals. This can be useful if the
+#' instrument is not available for the full sample period. Therefore only
+#' a part of the residuals corresponding to the length of instrument is used.
+#' For correct selection of residuals user is expected to specify rownames(data)
+#' and names(instrument) in same format. The subset is then based on string
+#' matching.
 #'
 #' @param horizon Integer scalar. The horizon for which impulse responses
 #' (and FEVDs) should be computed. Note that the first period corresponds to
@@ -26,7 +36,7 @@
 #' @param identification Logical scalar. Whether or not the shocks used for
 #' calculating impulses should be identified. Defaults to \code{TRUE}, i.e.
 #' identification via Cholesky decomposition of the VCOV-matrix unless
-#' \emph{sign_restr} is provided.
+#' \emph{sign_restr} or \emph{instrument} is provided.
 #' @param sign_restr Elements inform about expected impacts
 #' of certain shocks. Can be either \eqn{1}, \eqn{-1} or \eqn{0} depending
 #' on whether a positive, a negative or no contemporaneous effect of a
@@ -37,20 +47,45 @@
 #' \eqn{M - j} zero restrictions can be imposed on the \eqn{j}'th column.
 #' @param sign_lim Integer scalar. Maximum number of tries to find suitable
 #' matrices to for fitting sign or zero and sign restrictions.
+#' @param instrument Numeric vector. If provided, the identification is performed
+#' using proxy SVAR. The implemented proxy SVAR allows for different length of
+#' instrument and residuals. In such case a common subset is used for the
+#' identification. To subset correct observations names(instrument) and
+#'  rownames(data) are used - they need to have same format. See details
+#' @param proxyvar character string. Variable for which instrument is provided.
+#' Needs to correspond to one of the column names of the data used in
+#' \emph{bvar}.
+#' @param manual_matching If set to TRUE, no indices for the data and instrument
+#' are needed. Instead, the length of the instrument needs to
+#' match the length of the residuals. This rules out the possibility of performing
+#' identification on a proper subset of residuals.
 #'
 #' @return Returns a named list of class \code{bv_irf} with options for
 #' \code{\link{bvar}}, \code{\link{irf.bvar}} or \code{\link{fevd.bvar}}.
 #'
 #' @references
-#'   Rubio-Ramirez, J. F. and Waggoner, D. F. and Zha, T. (2010) Structural
-#'   Vector Autoregressions: Theory of Identification and Algorithms for
-#'   Inference. \emph{The Review of Economic Studies}, \bold{77}, 665-696,
-#'   \doi{10.1111/j.1467-937X.2009.00578.x}.
-#'   Arias, J.E. and Rubio-Ramirez, J. F. and Waggoner, D. F. (2018)
+#'  Arias, J.E. and Rubio-Ramirez, J. F. and Waggoner, D. F. (2018)
 #'   Inference Based on Structural Vector Autoregressions Identifiied with
 #'   Sign and Zero Restrictions: Theory and Applications.
 #'   \emph{Econometrica}, \bold{86}, 2, 685-720,
-#'   \doi{10.3982/ECTA14468}.
+#'   doi:10.3982/ECTA14468. \cr
+#'  Degasperi, R. and Ricco, G. (2021)
+#'   Information and policy shocks in monetary surprises. Working paper, University of Warwick. \cr
+#'  Mertens, K., and Ravn, M. O. (2013) The dynamic effects of personal and
+#'   corporate income tax changes in the United States. \emph{American
+#'   economic review}, \bold{103(4)}, 1212-1247, doi:10.1257/aer.103.4.1212. \cr
+#'  Miranda-Agrippino, S., & Ricco, G. (2021) The transmission of monetary
+#'   policy shocks. \emph{American Economic Journal: Macroeconomics},
+#'   \bold{13(3)}, 74-107, doi:10.1257/mac.20180124 \cr
+#'  Rubio-Ramirez, J. F. and Waggoner, D. F. and Zha, T. (2010) Structural
+#'   Vector Autoregressions: Theory of Identification and Algorithms for
+#'   Inference. \emph{The Review of Economic Studies}, \bold{77}, 665-696,
+#'   doi:10.1111/j.1467-937X.2009.00578.x. \cr
+#'  Stock, J. H., and Watson, M. W. (2012) Disentangling the Channels of the
+#'   2007-2009 Recession (No. w18094) \emph{National Bureau of Economic
+#'   Research}, doi:10.3386/w18094.
+
+
 #'
 #' @seealso \code{\link{irf.bvar}}; \code{\link{plot.bvar_irf}}
 #'
@@ -71,14 +106,32 @@
 #' zero_signs <- matrix(c(1, 0, NA, -1, 1, 0, -1, 1, 1), nrow = 3)
 #' bv_irf(sign_restr = zero_signs)
 #'
+#' # Set up structural impulse responses using external instrument with length
+#' corresponding to reduced form residuals
+#' bv_irf(instrument = istrument)
+#'
 #' # Prepare to estimate unidentified impulse responses
 #' bv_irf(identification = FALSE)
+#'
+#' # Prepare to estimate SVAR-IV using an indexed instrument for monetary policy
+#' # shocks. names(instrument) need to be same format as rownames(data), since
+#' # subsetting is based on string matching.
+#' data("mpi")
+#' bv_irf(instrument = instrument, manual_matching = FALSE, proxyvar = "FEDFUNDS")
+#'
+
+
+
 bv_irf <- function(
   horizon = 12,
   fevd = FALSE,
   identification = TRUE,
   sign_restr = NULL,
-  sign_lim = 1000) {
+  sign_lim = 1000,
+  instrument = NULL,
+  proxyvar = NULL,
+  manual_matching = FALSE
+  )  {
 
   # Input checks
   horizon <- int_check(horizon, min = 1, max = 1e6,
@@ -115,13 +168,43 @@ bv_irf <- function(
       }
     }
     # Cholesky
+
+  }
+
+  # IV
+
+  if(!is.null(instrument)){
+    if (is.numeric(instrument)) {
+    } else {
+      stop("Input must be a numeric vector. Multiple instruments are not
+           supported.")
+    }
+  }
+
+  if(!manual_matching && !is.null(instrument) &&
+     is.null(names(instrument))){
+    stop("No index detected for the object instrument. If manual_matching is
+  FALSE, names(instrument) need to be specified. For identification to work
+  correctly a common index for the data and the instrument must be specified.
+  Alternatively, switch manual_matching to TRUE.")
+  }
+
+  if (!is.null(instrument) && is.null(proxyvar)) {
+    stop("If 'instrument' is specified, 'proxyvar' must also be specified.")
+  }
+
+  if(!is.null(proxyvar) && !is.character(proxyvar)){
+    stop("Proxyvar needs to be an object of type character.")
   }
 
   # Outputs
   out <- list("horizon" = horizon, "fevd" = fevd,
     "identification" = identification,
     "sign_restr" = sign_restr, "zero" = zero,
-    "sign_lim" = sign_lim
+    "sign_lim" = sign_lim,
+    "instrument" = instrument,
+    "manual_matching" = manual_matching,
+    "proxyvar" = proxyvar
   )
 
   class(out) <- "bv_irf"
diff --git a/R/61_irf_compute.R b/R/61_irf_compute.R
index 46ded31..4b0e1f3 100644
--- a/R/61_irf_compute.R
+++ b/R/61_irf_compute.R
@@ -39,16 +39,50 @@ compute_irf <- function(
   M, lags,
   horizon,
   identification,
-  sign_restr, zero = FALSE, sign_lim = 10000) {
+  sign_restr, zero = FALSE, sign_lim = 10000,
+  residuals = NULL, instrument = NULL, manual_matching = FALSE, proxyvar = NULL) {
+
+
 
   # Identification
   if(identification) {
     sigma_chol <- t(chol(sigma))
-    if(is.null(sign_restr)) {
+    if(is.null(sign_restr) & is.null(instrument)) {
       shock <- sigma_chol
-    } else {
-      shock <- sign_restr(sigma_chol = sigma_chol,
-        sign_restr = sign_restr, M = M, sign_lim = sign_lim, zero = zero)
+    }
+    if(!is.null(sign_restr) & is.null(instrument)){
+    shock <- sign_restr(sigma_chol = sigma_chol,
+    sign_restr = sign_restr, M = M, sign_lim = sign_lim, zero = zero)
+    }
+    if(is.null(sign_restr) & !is.null(instrument)){
+          shock <- diag(M)
+
+          # Find the column index of the proxy variable
+          col_index <- which(colnames(residuals) == proxyvar)
+
+          # If the proxy variable is not already the first column, reorder columns
+          if (col_index != 1) {
+            # Instrumented variable swaps place with first variable in system
+            proxy_svar_ordering <- get_swapped_index(residuals, proxyvar)
+
+            # For proxy svar the instrumented variable is now ordered first
+            proxy_svar_output <- proxy_svar(
+              residuals[, proxy_svar_ordering],instrument)
+
+            # Reverse the ordering for output
+            proxy_svar_output$impact <- proxy_svar_output$impact[proxy_svar_ordering, , drop = FALSE]
+            shock[,col_index] <- proxy_svar_output$impact
+
+            } else{
+
+
+          # Assumed instrumented variable is ordered first
+          proxy_svar_output <- proxy_svar(residuals,instrument)
+          shock[,1] <- proxy_svar_output$impact
+            }
+    }
+    if(!is.null(sign_restr) & !is.null(instrument)){
+    stop("Sign restrictions and instrument cannot be used at the same time!")
     }
   } else {shock <- sigma}
 
@@ -60,5 +94,12 @@ compute_irf <- function(
   }
   irf_comp <- irf_comp[1:M, , 1:M]
 
-  return(irf_comp)
+  output <- list()
+  output$irf_comp <- irf_comp
+
+
+  output$iv_f_stat <- if(!is.null(instrument)){
+    proxy_svar_output$f_stat} else{NULL}
+
+  return(output)
 }
diff --git a/R/62b_proxy_var.R b/R/62b_proxy_var.R
new file mode 100644
index 0000000..a7a2557
--- /dev/null
+++ b/R/62b_proxy_var.R
@@ -0,0 +1,144 @@
+#' Proxy SVAR
+#'
+#' Performs identification using external \emph{instrument}. At this stage the
+#' \emph{instrument} and \emph{residuals} need to have equal number of
+#' observations. Transformations such as intersecting the instrument and
+#' residuals need to be done before calling this function.
+#'
+#' @param residuals Numeric matrix. During estimation, draw specific residuals
+#' are used.
+#' @param instrument Numeric vector containing the instrument. Possibly
+#' shortened beforehand to match the length of residuals.
+#'
+#' @return Returns a list with the following elements:
+#' \itemize{
+#'   \item \code{impact} - Contemporaneous transmission coefficients: Bzero
+#'   \item \code{Gamma}  - Estimated correlation between shock and instrument
+#'   \item \code{f_stat} - First stage F statistic
+#'   \item \code{lambda} - Reliability of instrument
+#' }
+#' @noRd
+
+proxy_svar <- function(residuals, instrument){
+
+  residuals <- as.matrix(residuals)
+
+  t <- nrow(residuals)
+  n <- ncol(residuals)
+  m <- 1 #only m=1 instrument implemented at the moment
+
+  # Create the Kronecker product
+  eye_n <- diag(rep(1, n)) #Matrix(diag(rep(1, n)), sparse = TRUE)
+  tempX <- cbind(1, instrument)
+  kron_matrix <- kronecker(eye_n, tempX)
+
+  # Solve the linear equation
+  betaIV <- solve(t(kron_matrix) %*% kron_matrix) %*% t(kron_matrix) %*% as.vector(residuals)
+
+  # Reshape betaIV
+  betaIV <- matrix(betaIV, nrow = length(betaIV) / n, byrow = FALSE)
+
+  # Transpose to get the final result
+  betaIV <- as.matrix(t(betaIV))
+
+  # F stat (regression on instruments of relevant innovations)
+  proxyVar <- instrument #Using notation from Agrip Ricco for this
+  tempU <- residuals[, 1:m] - tempX %*% betaIV[1:m, ]
+
+  tempY <- tempX %*% betaIV[1:m, ] - matrix(rep(mean(residuals[, 1:m]), t), ncol = m, byrow = TRUE)
+  k <- length(betaIV[1:m, ]) - 1
+
+  f_stat <- ((t(tempY) %*% tempY) / k) / ((t(tempU) %*% tempU) / (t - k - 1))
+
+
+  ################
+  # Assuming m is defined and betaIV is already a numeric matrix [1:5, 1:2]
+
+  # beta_11 and beta_21
+  beta_11 <- betaIV[1:m, 2:(m+1)]
+  beta_21 <- betaIV[(m + 1):nrow(betaIV), 2:(m+1)]
+
+  # ratio of regression coefficients
+  B21B11 <- beta_21 / beta_11
+
+  # Covariance matrix
+  SigmaU <- cov(residuals)
+
+  # Identification
+  Zeta <- (B21B11 * SigmaU[1:m, 1:m]) %*% t(B21B11) -
+    (SigmaU[(m + 1):nrow(SigmaU), 1:m] %*% t(B21B11) + B21B11 %*% t(SigmaU[(m + 1):nrow(SigmaU), 1:m])) +
+    SigmaU[(m + 1):nrow(SigmaU), (m + 1):ncol(SigmaU)]
+
+  B12B12 <- t(SigmaU[(m + 1):nrow(SigmaU), 1:m] - B21B11 * SigmaU[1:m, 1:m]) %*% solve(Zeta) %*%
+    (SigmaU[(m + 1):nrow(SigmaU), 1:m] - B21B11 * SigmaU[1:m, 1:m])
+
+  B11B11 <- SigmaU[1:m, 1:m] - B12B12
+
+  #Since we only have 1 instrument
+  B11 <- sqrt(B11B11)   # beta_{11}
+  B <- as.vector(B11) * c(1, B21B11)  # first column of B (u_t = B * e_t)
+
+  ###########################
+  #%realized shock sequences (Montiel-Olea, Stock and Watson)
+  # Assuming T, residuals, and proxyVar are defined
+
+  # Create tempX matrix
+  tempX <- cbind(1, residuals)
+
+  # Calculate e using matrix multiplication and solving the system of equations
+  e <-  tempX %*% solve(t(tempX) %*% tempX) %*% t(tempX) %*% proxyVar
+
+  # Standardize e to have unit variance
+  e <- scale(e, center = TRUE, scale = apply(e, 2, sd)) # Now, e contains the unit variance shock series
+
+  ##################
+  #Again since m=1
+
+  # Assuming proxyVar, residuals, and e are defined in R
+
+  # Proportion of uncensored data
+  D <- (proxyVar != 0)
+
+  # Covariance matrix
+  SigmaMU <- cov(cbind(proxyVar[D], residuals[D, ]), use = "pairwise")
+
+  # Relevance
+  A_ <- SigmaMU[1,2:ncol(SigmaMU)]
+  Phi <- 1/(MASS::ginv(A_) %*% B)
+  G <- 1/(t %*% MASS::ginv(sum(D)))
+  Gamma <- (Phi*(1/G))
+
+  # Calculate eSquare and zSquare
+  eSquare <- e^2
+  zSquare <- (proxyVar - rep(Gamma,t) * as.vector(e))^2
+
+  # Calculate Lambda
+  Lambda <- 1/((Gamma^2 * sum(eSquare[D]) + sum(zSquare[D])) / (Gamma^2 * sum(eSquare[D])))
+
+
+  ####################
+  iP <- 1 #instrument position
+  # Load output structure
+
+  impact = matrix(NA, n, m)      # contemporaneous transmission coefficients: Bzero
+  impact[iP,] <- B[1:m]
+  impact[(iP+1):n,] <- B[(m + 1):n]
+
+  #normalize
+  impact <- impact / impact[1]
+
+
+  #Summary
+  # impact = impact,             # contemporaneous transmission coefficients: Bzero
+  # Gamma = Gamma,             # estimated correlation between shock and instrument
+  # L = Lambda,                # reliability of instrument
+  # e = e,                     # realized shocks series
+  # fstat = diag(f_stat)       # F statistic of regression on instrument
+
+  output <- list(impact, Gamma, f_stat, Lambda)
+  names(output) <- c("impact", "Gamma", "f_stat", "Lambda")
+
+  return(output)
+}
+
+
diff --git a/R/64_irf_method.R b/R/64_irf_method.R
index 093e649..75ff8c3 100644
--- a/R/64_irf_method.R
+++ b/R/64_irf_method.R
@@ -39,13 +39,13 @@
 #'
 #' @examples
 #' \donttest{
-#' # Access a subset of the fred_qd dataset
-#' data <- fred_qd[, c("CPIAUCSL", "UNRATE", "FEDFUNDS")]
-#' # Transform it to be stationary
-#' data <- fred_transform(data, codes = c(5, 5, 1), lag = 4)
+#' # Access the fred_md dataset
+#' data <- fred_md[,c("CPIAUCSL", "UNRATE", "FEDFUNDS")]
+#' # Transform and subset
+#' data <- fred_transform(data, codes = c(4, 1, 1))[1:732,]
 #'
-#' # Estimate a BVAR using one lag, default settings and very few draws
-#' x <- bvar(data, lags = 1, n_draw = 600L, n_burn = 100L, verbose = FALSE)
+#' # Estimate a BVAR using 12 lags, default settings and very few draws
+#' x <- bvar(data, lags = 12, n_draw = 600L, n_burn = 100L, verbose = FALSE)
 #'
 #' # Compute + store IRF with a longer horizon, no identification and thinning
 #' irf(x) <- irf(x, bv_irf(horizon = 24L, identification = FALSE), n_thin = 5L)
@@ -67,6 +67,11 @@
 #'
 #' # Limit the summary to responses of variable #2
 #' summary(x, vars_response = 2L)
+#'
+#' # Compute IRF using an external instrument for monetary policy shocks
+#' data("mpi")
+#' irf(x) <- irf.bvar(x, bv_irf(horizon = 24L, instrument = instrument,
+#'                              proxyvar = "FEDFUNDS"), n_thin = 1L)
 #' }
 irf.bvar <- function(x, ..., conf_bands, n_thin = 1L, verbose = FALSE) {
 
@@ -74,7 +79,6 @@ irf.bvar <- function(x, ..., conf_bands, n_thin = 1L, verbose = FALSE) {
   irf_store <- x[["irf"]]
   verbose <- isTRUE(verbose)
 
-
   # Calculate impulse responses -----
 
   if(is.null(irf_store) || length(dots) != 0L) {
@@ -84,12 +88,16 @@ irf.bvar <- function(x, ..., conf_bands, n_thin = 1L, verbose = FALSE) {
     irf <- if(length(dots) > 0 && inherits(dots[[1]], "bv_irf")) {
       dots[[1]]
     } else {bv_irf(...)}
+    instrument <- irf[["instrument"]]
+    manual_matching <- irf[["manual_matching"]]
+    proxyvar <- irf[["proxyvar"]]
 
     n_pres <- x[["meta"]][["n_save"]]
     n_thin <- int_check(n_thin, min = 1, max = (n_pres / 10),
       "Issue with n_thin. Maximum allowed is n_save / 10.")
     n_save <- int_check((n_pres / n_thin), min = 1)
 
+    X <- x[["meta"]][["X"]]
     Y <- x[["meta"]][["Y"]]
     N <- x[["meta"]][["N"]]
     K <- x[["meta"]][["K"]]
@@ -97,6 +105,11 @@ irf.bvar <- function(x, ..., conf_bands, n_thin = 1L, verbose = FALSE) {
     lags <- x[["meta"]][["lags"]]
     beta <- x[["beta"]]
     sigma <- x[["sigma"]]
+    #For the case we need to reorder due to using IV
+    vars <- name_deps(x[["variables"]], M = M)
+    vars_expl <- name_expl(vars, M = M, lags = lags)
+    dimnames(beta) <- list(NULL, vars_expl, vars)
+
 
     # Check sign restrictions
     if(!is.null(irf[["sign_restr"]]) && length(irf[["sign_restr"]]) != M ^ 2 ||
@@ -104,6 +117,8 @@ irf.bvar <- function(x, ..., conf_bands, n_thin = 1L, verbose = FALSE) {
       stop("Dimensions of provided restrictions do not fit the data.")
     }
 
+
+
     # Sampling ---
 
     irf_store <- structure(list(
@@ -112,7 +127,10 @@ irf.bvar <- function(x, ..., conf_bands, n_thin = 1L, verbose = FALSE) {
         structure(
           list("fevd" = array(NA, c(n_save, M, irf[["horizon"]], M)),
             "variables" = x[["variables"]]), class = "bvar_fevd")
-      } else {NULL}, "setup" = irf, "variables" = x[["variables"]]),
+      } else {NULL}, "setup" = irf, "variables" = x[["variables"]],
+      "iv_stats" = if(!is.null(instrument)){
+        list("f_stat_store" = array(NA, n_save), "f_stat" = list())
+        } else {NULL}),
       class = "bvar_irf")
 
     j <- 1
@@ -122,12 +140,32 @@ irf.bvar <- function(x, ..., conf_bands, n_thin = 1L, verbose = FALSE) {
     }
     for(i in seq_len(n_save)) {
       beta_comp <- get_beta_comp(beta[j, , ], K, M, lags)
-      irf_comp  <- compute_irf(
+
+      if(!is.null(instrument)){
+        # For identification, if IV is shorter than residuals, subset residuals.
+        #From 62b_proxy_var.R
+        residuals_draw <- Y - X %*% beta[j, , ]
+        intersection <- check_iv(residuals_draw, instrument,
+                                 manual_matching = manual_matching,
+                                 proxyvar = proxyvar)
+        if(i==1){print(intersection)}
+      } else{
+        intersection <- list()
+        intersection$residuals_shortened <- NULL
+        intersection$instrument_shortened <- NULL
+      }
+      output  <- compute_irf(
         beta_comp = beta_comp, sigma = sigma[j, , ], M = M, lags = lags,
         horizon = irf[["horizon"]], identification = irf[["identification"]],
         sign_restr = irf[["sign_restr"]], zero = irf[["zero"]],
-        sign_lim = irf[["sign_lim"]])
+        sign_lim = irf[["sign_lim"]],
+        residuals = intersection$residuals_shortened,
+        instrument = intersection$instrument_shortened,
+        manual_matching = manual_matching,
+        proxyvar = proxyvar)
+      irf_comp <- output$irf_comp
       irf_store[["irf"]][i, , , ] <- irf_comp
+      irf_store[["iv_stats"]][["f_stat_store"]][i] <- output$iv_f_stat
 
       if(irf[["fevd"]]) { # Forecast error variance decomposition
         irf_store[["fevd"]][["fevd"]][i, , , ] <- compute_fevd(
@@ -158,6 +196,13 @@ irf.bvar <- function(x, ..., conf_bands, n_thin = 1L, verbose = FALSE) {
   }
 
 
+  f_stats <- quantile(irf_store[["iv_stats"]][["f_stat_store"]], c(0.025, 0.05,0.5,0.95, 0.975))
+
+  irf_store[["iv_stats"]][["f_stat"]] <- list(
+    "quantiles" = round(f_stats,2)
+    #"conf_interval (90%)" = c(lower = 1.8, upper = 3.2)
+  )
+
   return(irf_store)
 }
 
diff --git a/R/65_irf_print.R b/R/65_irf_print.R
index e3c6ef3..885b4a8 100644
--- a/R/65_irf_print.R
+++ b/R/65_irf_print.R
@@ -51,8 +51,19 @@ print.bvar_fevd <- function(x, digits = 4L, complete = FALSE, ...) {
   cat("\nIdentification: ")
   if(x[["identification"]]) {
     if(is.null(x[["sign_restr"]])) {
+      if(is.null(x[["instrument"]])) {
       cat("Cholesky decomposition")
-    } else {
+      } else{cat("External Instrument"
+                 # ,
+                 # "\nFirst stage F-statistic (95% CI): ",
+                 # x[["irf"]][["iv_stats"]][["f_stat"]]$quantiles[3],
+                 # "(",
+                 # x[["irf"]][["iv_stats"]][["f_stat"]]$quantiles[1], ",",
+                 # x[["irf"]][["iv_stats"]][["f_stat"]]$quantiles[5],
+                 # ")"
+                 )
+        }
+    } else{
       cat("Sign restrictions", "\nChosen restrictions:\n", sep = "")
       sign_restr <- apply(x[["sign_restr"]], 2, factor,
         levels = c(-1, 0, 1), labels = c("-", "0", "+"))
diff --git a/R/90_print.R b/R/90_print.R
index 843306b..06a5058 100644
--- a/R/90_print.R
+++ b/R/90_print.R
@@ -19,3 +19,46 @@ print.bvar <- function(x, ...) {
 
   return(invisible(x))
 }
+
+#' @export
+print.check_iv <- function(out, ...) {
+
+  # Extract the start and end dates from the data and instrument date ranges
+  start_date_data <- out$dates_residuals[1]
+  end_date_data <- tail(out$dates_residuals, 1)
+  start_date_instrument <- out$dates_instrument[1]
+  end_date_instrument <- tail(out$dates_instrument, 1)
+
+  if(out$manual_matching == TRUE){
+    #print(class(out$residuals_shortened))
+    cat("Found residuals of length", nrow(out$residuals_shortened),  "\n")
+    cat("Found instrument of length", length(out$instrument_shortened),  "\n")
+    # cat("Head of cbind(residuals[,1], instrument):", "\n")
+    # m1 <- cbind(out$residuals_shortened[,1], out$instrument_shortened)
+    # display_matrix(m1)
+    #print(head(cbind(out$residuals_shortened[,1], out$instrument_shortened), 3))
+    #print(tail(cbind(out$residuals_shortened[,1], out$instrument_shortened), 3))
+    # cat("Found residuals:", head(irf$residuals_shortened,1), "...", tail(out$residuals_shortened,1),
+    #     "and instrument:",  head(out$instrument_shortened,1), "...",  tail(out$instrument_shortened,1))
+  } else{
+
+  # Get the number of matched dates and display the first and last matched observations
+  num_matches <- length(out$common_dates)
+  matched_display <- paste0(
+    out$common_dates[1], " (instrument value ", format(out$instrument_shortened[out$common_dates[1]], scientific = TRUE),
+    ") until ",
+    out$common_dates[num_matches], " (", format(out$instrument_shortened[out$common_dates[num_matches]], scientific = TRUE),
+    ")"
+  )
+
+  # Print the summary information
+  cat("Found residuals from", start_date_data, "until", end_date_data,
+      "and instrument from", start_date_instrument, "until", end_date_instrument, "\n")
+  cat("Matched", num_matches, "observations: from ", matched_display, "\n")
+
+   }
+
+  return(invisible(out))
+}
+
+
diff --git a/R/data_mpi.R b/R/data_mpi.R
new file mode 100644
index 0000000..0a7e92c
--- /dev/null
+++ b/R/data_mpi.R
@@ -0,0 +1,28 @@
+
+#' Monetary Policy Instrument
+#'
+#' We make use of time series at monthly frequency of the high-frequency
+#' instruments for monetary policy and information shocks as defined in
+#' Miranda-Agrippino and Ricco (2021) and extended to 2018:12 in
+#' Degasperi and Ricco (2021).
+
+#' The monetary policy instrument that is provided here is licensed
+#' under a modified CC-BY 4.0 license that can be found in the provided
+#' \emph{LICENSE} file. #' For further details see
+#' \url{https://github.com/riccardo-degasperi/info-policy-surprises}.
+#' We would like to thank Silvia Miranda-Agrippino, Giovanni Ricco and
+#' Riccardo Degasperi for creating, updating and making the instrument available.
+#'
+#' @name mpi
+#' @docType data
+#' @usage data(mpi)
+#' @format A \code{numeric vector} object with dates as rownames.
+#' @keywords monetary policy; shocks; instrument
+#' @references
+#'   Miranda-Agrippino, S., & Ricco, G. (2021). The transmission of monetary
+#'   policy shocks. \emph{American Economic Journal: Macroeconomics},
+#'    \bold{13(3)}, 74-107, doi:10.1257/mac.20180124 \cr
+#'    Degasperi, R. and Ricco, G., 2021.
+#'    Information and policy shocks in monetary surprises. Working Paper.
+#' @source \url{https://github.com/riccardo-degasperi/info-policy-surprises}
+"mpi"
diff --git a/data/fred_md.rda b/data/fred_md.rda
index fb7d9fc..7c38ffb 100644
Binary files a/data/fred_md.rda and b/data/fred_md.rda differ
diff --git a/data/fred_permitted.txt b/data/fred_permitted.txt
index 1a3cda9..59b82de 100644
--- a/data/fred_permitted.txt
+++ b/data/fred_permitted.txt
@@ -153,10 +153,8 @@ GS1TB3Mx
 GS10TB3Mx
 CPF3MTB3Mx
 BOGMBASEREALx
-IMFSLx
 M1REAL
 M2REAL
-MZMREAL
 BUSLOANSx
 CONSUMERx
 NONREVSLx
@@ -172,9 +170,9 @@ NWPIx
 TARESAx
 HNOREMQ027Sx
 TFAABSHNOx
-VXOCLSx
+VIXCLSx
 USSTHPI
-TWEXMMTH
+TWEXAFEGSMTHx
 EXUSEU
 EXSZUSx
 EXJPUSx
@@ -235,6 +233,9 @@ TABSNNBx
 TNWBSNNBx
 TNWBSNNBBDIx
 CNCFx
+S.P.500
+S.P.div.yield
+S.P.PE.ratio
 RPI
 W875RX1
 DPCERA3M086SBEA
@@ -264,4 +265,3 @@ DNDGRG3M086SBEA
 DSERRG3M086SBEA
 CES2000000008
 CES3000000008
-MZMSL
diff --git a/data/fred_qd.rda b/data/fred_qd.rda
index 7350cee..cc2dc2b 100644
Binary files a/data/fred_qd.rda and b/data/fred_qd.rda differ
diff --git a/data/mpi.rda b/data/mpi.rda
new file mode 100644
index 0000000..adace1a
Binary files /dev/null and b/data/mpi.rda differ
diff --git a/inst/fred_trans.csv b/inst/fred_trans.csv
index dc55f01..0c32ecb 100644
--- a/inst/fred_trans.csv
+++ b/inst/fred_trans.csv
@@ -73,9 +73,8 @@
 "NONREVSL","log-2nd-diff",
 "CONSPI","1st-diff",
 "SP500","log-diff",
-"SPINDUST","log-diff",
-"SPDIVYIELD","1st-diff",
-"SPPERATIO","log-diff",
+"SPINDUST","1st-diff",
+"SPDIVYIELD","log-diff",
 "FEDFUNDS","1st-diff","1st-diff"
 "CP3Mx","1st-diff",
 "TB3MS","1st-diff","1st-diff"
@@ -276,6 +275,5 @@
 "TNWBSNNBBDIx",,"1st-diff"
 "CNCFx",,"log-diff"
 "S.P.500",,"log-diff"
-"S.P..indust",,"log-diff"
 "S.P.div.yield",,"1st-diff"
 "S.P.PE.ratio",,"log-diff"
diff --git a/inst/tinytest/test_BVAR.R b/inst/tinytest/test_BVAR.R
index 831ede4..8e7e2f8 100644
--- a/inst/tinytest/test_BVAR.R
+++ b/inst/tinytest/test_BVAR.R
@@ -1,7 +1,15 @@
-
+devtools::load_all()
 # API tests -------
 
-data <- data2 <- data3 <- matrix(rnorm(1000), nrow = 200)
+data <- data2 <- data3 <- data_iv <- matrix(rnorm(1000), nrow = 200)
+
+colnames(data_iv) <- paste0("col", 1:ncol(data_iv))
+proxyvar <- colnames(data_iv)[2]
+
+instrument <-  rnorm(200)
+instrument_named <- instrument
+names(instrument_named) <- seq(1:200)
+instrument_exact <- instrument[3:200]
 
 # Fail and prepare -----
 
@@ -120,11 +128,29 @@ expect_silent(opt_irf3 <- bv_irf(fevd = FALSE, # Zero sign restricted
 expect_silent(bv_irf(sign_restr = c(1, NA, -1, 1), sign_lim = 1000))
 expect_silent(bv_irf(sign_restr = c(0, NA, NA, 1), sign_lim = 1000))
 
-# Underidentified, too many 0, non-square restrictions
+expect_silent(opt_irf4_iv <- bv_irf(instrument = instrument,
+                                    manual_matching = TRUE,
+                                    proxyvar = proxyvar))
+expect_silent(opt_irf4_iv_named <- bv_irf(instrument = instrument_named,
+                                          proxyvar = proxyvar))
+expect_silent(opt_irf4_iv_exact <- bv_irf(instrument = instrument_exact,
+                                 manual_matching = TRUE,
+                                 proxyvar = proxyvar))
+expect_silent(print(opt_irf4_iv))
+expect_silent(print(opt_irf4_iv_named))
+expect_silent(print(opt_irf4_iv_exact))
+
+
+# Underidentified, too many 0, non-square restrictions, wrong input type
 expect_message(bv_irf(sign_restr = matrix(c(NA, NA, NA, NA), nrow = 2)))
 expect_error(bv_irf(sign_restr = matrix(c(0, 0, -1, NA), nrow = 2)))
 expect_error(bf_irf(sign_restr = matrix(rnorm(6), nrow = 3)))
 
+# wrong input type, missing index or missing mandatory inputs
+expect_error(bv_irf(instrument = data.frame(instrument)))
+expect_error(bv_irf(instrument = instrument))
+expect_error(bv_irf(instrument = instrument_named))
+
 
 # Run and analyse -----
 
@@ -139,6 +165,9 @@ expect_silent(run2 <- bvar(data[, 1:3], lags = 2,
 expect_silent(run3 <- bvar(data, lags = 2,
   priors = bv_priors(hyper = c("lambda", "psi")), mh = mh,
   n_draw = 1000L, n_burn = 500L))
+expect_silent(run4 <- bvar(data_iv, lags = 2,
+                      priors = bv_priors(hyper = c("lambda", "psi")), mh = mh,
+                      n_draw = 1000L, n_burn = 500L))
 
 
 # 5*_fcast ---
@@ -161,6 +190,8 @@ expect_silent(irf(run) <- irf(run, opt_irf1))
 expect_silent(irfs1 <- irf(run, verbose = TRUE))
 expect_silent(irfs2 <- irf(run2, opt_irf2))
 expect_silent(irfs3 <- irf(run2, opt_irf3))
+expect_silent(irfs4_iv_exact <- irf(run4, opt_irf4_iv_exact))
+expect_silent(irfs4_iv_named <- irf(run4, opt_irf4_iv_named))
 
 expect_silent(print(irfs1))
 expect_silent(print(summary(irfs1)))
@@ -169,6 +200,13 @@ expect_silent(print(fevd(run2))) # Recalculates
 expect_silent(print(fevd(irfs2))) # Recalculates
 expect_silent(plot(irfs1, vars_res = 1, vars_imp = 1))
 
+# No names / wrong names provided and exact_matching = F
+expect_error(irf(run4, bv_irf(instrument = instrument,
+                              manual_matching = F,
+                              proxyvar = proxyvar)))
+# manual_matching with wrong length
+expect_error(irf(run4, opt_irf4_iv))
+
 
 # 80_coda ---
 
diff --git a/man/BVAR-package.Rd b/man/BVAR-package.Rd
index fa7be2e..75aff67 100644
--- a/man/BVAR-package.Rd
+++ b/man/BVAR-package.Rd
@@ -43,6 +43,7 @@ Authors:
 Other contributors:
 \itemize{
   \item Nirai Tomass [contributor]
+  \item Gabriel Konecny [contributor]
   \item Michael McCracken [data contributor]
   \item Serena Ng [data contributor]
 }
diff --git a/man/bv_dummy.Rd b/man/bv_dummy.Rd
index 2421c65..1e2ea72 100644
--- a/man/bv_dummy.Rd
+++ b/man/bv_dummy.Rd
@@ -6,11 +6,11 @@
 \alias{bv_sur}
 \title{Dummy prior settings}
 \usage{
-bv_dummy(mode = 1, sd = 1, min = 0.0001, max = 5, fun)
+bv_dummy(mode = 1, sd = 1, min = 1e-04, max = 5, fun)
 
-bv_soc(mode = 1, sd = 1, min = 0.0001, max = 50)
+bv_soc(mode = 1, sd = 1, min = 1e-04, max = 50)
 
-bv_sur(mode = 1, sd = 1, min = 0.0001, max = 50)
+bv_sur(mode = 1, sd = 1, min = 1e-04, max = 50)
 }
 \arguments{
 \item{mode, sd}{Numeric scalar. Mode / standard deviation of the
diff --git a/man/bv_irf.Rd b/man/bv_irf.Rd
index 59d938d..f41aa72 100644
--- a/man/bv_irf.Rd
+++ b/man/bv_irf.Rd
@@ -9,7 +9,10 @@ bv_irf(
   fevd = FALSE,
   identification = TRUE,
   sign_restr = NULL,
-  sign_lim = 1000
+  sign_lim = 1000,
+  instrument = NULL,
+  proxyvar = NULL,
+  manual_matching = FALSE
 )
 }
 \arguments{
@@ -23,7 +26,7 @@ decompositions should be calculated.}
 \item{identification}{Logical scalar. Whether or not the shocks used for
 calculating impulses should be identified. Defaults to \code{TRUE}, i.e.
 identification via Cholesky decomposition of the VCOV-matrix unless
-\emph{sign_restr} is provided.}
+\emph{sign_restr} or \emph{instrument} is provided.}
 
 \item{sign_restr}{Elements inform about expected impacts
 of certain shocks. Can be either \eqn{1}, \eqn{-1} or \eqn{0} depending
@@ -36,6 +39,21 @@ value is \code{NULL}. Note that in order to be fully identified at least
 
 \item{sign_lim}{Integer scalar. Maximum number of tries to find suitable
 matrices to for fitting sign or zero and sign restrictions.}
+
+\item{instrument}{Numeric vector. If provided, the identification is performed
+using proxy SVAR. The implemented proxy SVAR allows for different length of
+instrument and residuals. In such case a common subset is used for the
+identification. To subset correct observations names(instrument) and
+ rownames(data) are used - they need to have same format. See details}
+
+\item{proxyvar}{character string. Variable for which instrument is provided.
+Needs to correspond to one of the column names of the data used in
+\emph{bvar}.}
+
+\item{manual_matching}{If set to TRUE, no indices for the data and instrument
+are needed. Instead, the length of the instrument needs to
+match the length of the residuals. This rules out the possibility of performing
+identification on a proper subset of residuals.}
 }
 \value{
 Returns a named list of class \code{bv_irf} with options for
@@ -49,16 +67,26 @@ or not forecast error variance decompositions (FEVDs) should be included
 as well as if and what kind of identification should be used. See the Details
 section for further information on identification. Identification can be
 achieved via Cholesky decomposition, sign restrictions (Rubio-Ramirez,
-Waggoner and Zha, 2010), and zero and sign restrictions (Arias,
-Rubio-Ramirez and Waggoner, 2018).
+Waggoner and Zha, 2010), combination of zero and sign restrictions (Arias,
+Rubio-Ramirez and Waggoner, 2018) or using an external instrument
+(Stock and Watson 2012, Mertens and Ravn 2013)
 }
 \details{
 Identification can be performed via Cholesky decomposition, sign
-restrictions, or zero and sign restrictions. The algorithm
-for generating suitable sign restrictions follows Rubio-Ramirez, Waggoner
-and Zha (2010), while the one for zero and sign restrictions follows
-Arias, Rubio-Ramirez and Waggoner (2018).
+restrictions, zero and sign restrictions or using an external instrument.
+The algorithm for generating suitable sign restrictions follows
+Rubio-Ramirez, Waggoner and Zha (2010), while the one for zero and sign
+restrictions follows Arias, Rubio-Ramirez and Waggoner (2018).
 Note the possiblity of finding no suitable zero/sign restrictions.
+The proxy SVAR identification follows the implementation in
+Miranda-Agrippino and Ricco (2021) and examples make use of their instrument
+extended in Degasperi and Ricco (2021). The implemented proxy SVAR allows
+for identification using subsample of residuals. This can be useful if the
+instrument is not available for the full sample period. Therefore only
+a part of the residuals corresponding to the length of instrument is used.
+For correct selection of residuals user is expected to specify rownames(data)
+and names(instrument) in same format. The subset is then based on string
+matching.
 }
 \examples{
 # Set impulse responses to a horizon of 20 time periods and enable FEVD
@@ -73,19 +101,41 @@ bv_irf(sign_restr = signs)
 zero_signs <- matrix(c(1, 0, NA, -1, 1, 0, -1, 1, 1), nrow = 3)
 bv_irf(sign_restr = zero_signs)
 
+# Set up structural impulse responses using external instrument with length
+corresponding to reduced form residuals
+bv_irf(instrument = istrument)
+
 # Prepare to estimate unidentified impulse responses
 bv_irf(identification = FALSE)
+
+# Prepare to estimate SVAR-IV using an indexed instrument for monetary policy
+# shocks. names(instrument) need to be same format as rownames(data), since
+# subsetting is based on string matching.
+data("mpi")
+bv_irf(instrument = instrument, manual_matching = FALSE, proxyvar = "FEDFUNDS")
+
 }
 \references{
-Rubio-Ramirez, J. F. and Waggoner, D. F. and Zha, T. (2010) Structural
-  Vector Autoregressions: Theory of Identification and Algorithms for
-  Inference. \emph{The Review of Economic Studies}, \bold{77}, 665-696,
-  \doi{10.1111/j.1467-937X.2009.00578.x}.
-  Arias, J.E. and Rubio-Ramirez, J. F. and Waggoner, D. F. (2018)
+Arias, J.E. and Rubio-Ramirez, J. F. and Waggoner, D. F. (2018)
   Inference Based on Structural Vector Autoregressions Identifiied with
   Sign and Zero Restrictions: Theory and Applications.
   \emph{Econometrica}, \bold{86}, 2, 685-720,
-  \doi{10.3982/ECTA14468}.
+  doi:10.3982/ECTA14468. \cr
+ Degasperi, R. and Ricco, G., 2021.
+  Information and policy shocks in monetary surprises. Working Paper. \cr
+ Mertens, K., and Ravn, M. O. (2013). The dynamic effects of personal and
+  corporate income tax changes in the United States. \emph{American
+  economic review}, \bold{103(4)}, 1212-1247, doi:10.1257/aer.103.4.1212. \cr
+ Miranda-Agrippino, S., & Ricco, G. (2021). The transmission of monetary
+  policy shocks. \emph{American Economic Journal: Macroeconomics},
+  \bold{13(3)}, 74-107, doi:10.1257/mac.20180124 \cr
+ Rubio-Ramirez, J. F. and Waggoner, D. F. and Zha, T. (2010) Structural
+  Vector Autoregressions: Theory of Identification and Algorithms for
+  Inference. \emph{The Review of Economic Studies}, \bold{77}, 665-696,
+  doi:10.1111/j.1467-937X.2009.00578.x. \cr
+ Stock, J. H., and Watson, M. W. (2012). Disentangling the Channels of the
+  2007-2009 Recession (No. w18094) \emph{National Bureau of Economic
+  Research}, doi:10.3386/w18094.
 }
 \seealso{
 \code{\link{irf.bvar}}; \code{\link{plot.bvar_irf}}
diff --git a/man/bv_minnesota.Rd b/man/bv_minnesota.Rd
index dd32604..1222ff9 100644
--- a/man/bv_minnesota.Rd
+++ b/man/bv_minnesota.Rd
@@ -12,7 +12,7 @@ bv_minnesota(
   lambda = bv_lambda(),
   alpha = bv_alpha(),
   psi = bv_psi(),
-  var = 10000000,
+  var = 1e+07,
   b = 1
 )
 
@@ -20,11 +20,11 @@ bv_mn(
   lambda = bv_lambda(),
   alpha = bv_alpha(),
   psi = bv_psi(),
-  var = 10000000,
+  var = 1e+07,
   b = 1
 )
 
-bv_lambda(mode = 0.2, sd = 0.4, min = 0.0001, max = 5)
+bv_lambda(mode = 0.2, sd = 0.4, min = 1e-04, max = 5)
 
 bv_alpha(mode = 2, sd = 0.25, min = 1, max = 3)
 
diff --git a/man/fred_qd.Rd b/man/fred_qd.Rd
index 9a08eea..1233375 100644
--- a/man/fred_qd.Rd
+++ b/man/fred_qd.Rd
@@ -8,7 +8,7 @@
 \format{
 A \code{data.frame} object with dates as rownames.
 
-An object of class \code{data.frame} with 777 rows and 118 columns.
+An object of class \code{data.frame} with 776 rows and 126 columns.
 }
 \source{
 \url{https://research.stlouisfed.org/econ/mccracken/fred-databases/}
diff --git a/man/irf.bvar.Rd b/man/irf.bvar.Rd
index 8fca29f..5d065d8 100644
--- a/man/irf.bvar.Rd
+++ b/man/irf.bvar.Rd
@@ -73,13 +73,13 @@ confidence bands.
 }
 \examples{
 \donttest{
-# Access a subset of the fred_qd dataset
-data <- fred_qd[, c("CPIAUCSL", "UNRATE", "FEDFUNDS")]
-# Transform it to be stationary
-data <- fred_transform(data, codes = c(5, 5, 1), lag = 4)
+# Access the fred_md dataset
+data <- fred_md[,c("CPIAUCSL", "UNRATE", "FEDFUNDS")]
+# Transform and subset
+data <- fred_transform(data, codes = c(4, 1, 1))[1:732,]
 
-# Estimate a BVAR using one lag, default settings and very few draws
-x <- bvar(data, lags = 1, n_draw = 600L, n_burn = 100L, verbose = FALSE)
+# Estimate a BVAR using 12 lags, default settings and very few draws
+x <- bvar(data, lags = 12, n_draw = 600L, n_burn = 100L, verbose = FALSE)
 
 # Compute + store IRF with a longer horizon, no identification and thinning
 irf(x) <- irf(x, bv_irf(horizon = 24L, identification = FALSE), n_thin = 5L)
@@ -101,6 +101,11 @@ summary(x)
 
 # Limit the summary to responses of variable #2
 summary(x, vars_response = 2L)
+
+# Compute IRF using an external instrument for monetary policy shocks
+data("mpi")
+irf(x) <- irf.bvar(x, bv_irf(horizon = 24L, instrument = instrument,
+                             proxyvar = "FEDFUNDS"), n_thin = 1L)
 }
 }
 \seealso{
diff --git a/man/mpi.Rd b/man/mpi.Rd
new file mode 100644
index 0000000..14ed8e7
--- /dev/null
+++ b/man/mpi.Rd
@@ -0,0 +1,38 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/data_mpi.R
+\docType{data}
+\name{mpi}
+\alias{mpi}
+\title{Monetary Policy Instrument}
+\format{
+A \code{numeric vector} object with dates as rownames.
+}
+\source{
+\url{https://github.com/riccardo-degasperi/info-policy-surprises}
+}
+\usage{
+data(mpi)
+}
+\description{
+We make use of time series at monthly frequency of the high-frequency
+instruments for monetary policy and information shocks as defined in
+Miranda-Agrippino and Ricco (2021) and extended to 2018:12 in
+Degasperi and Ricco (2021).
+The monetary policy instrument that is provided here is licensed
+under a modified CC-BY 4.0 license that can be found in the provided
+\emph{LICENSE} file. #' For further details see
+\url{https://github.com/riccardo-degasperi/info-policy-surprises}.
+We would like to thank Silvia Miranda-Agrippino, Giovanni Ricco and
+Riccardo Degasperi for creating, updating and making the instrument available.
+}
+\references{
+Miranda-Agrippino, S., & Ricco, G. (2021). The transmission of monetary
+  policy shocks. \emph{American Economic Journal: Macroeconomics},
+   \bold{13(3)}, 74-107, doi:10.1257/mac.20180124 \cr
+   Degasperi, R. and Ricco, G., 2021.
+   Information and policy shocks in monetary surprises. Working Paper.
+}
+\keyword{instrument}
+\keyword{monetary}
+\keyword{policy;}
+\keyword{shocks;}
diff --git a/man/set_dates.Rd b/man/set_dates.Rd
new file mode 100644
index 0000000..252280a
--- /dev/null
+++ b/man/set_dates.Rd
@@ -0,0 +1,26 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/12_aux.R
+\name{set_dates}
+\alias{set_dates}
+\title{Set dates}
+\usage{
+set_dates(data, start, frequency)
+}
+\arguments{
+\item{data}{Numeric vector, matrix or a data.frame. Note that observations
+are expected to be ordered from earliest to latest, and variables in the
+columns.}
+
+\item{start}{Starting date of the data in form "YYYY-MM-DD".}
+
+\item{frequency}{Frequency of the data: "year", "quarter", "month", "week" or "day".}
+}
+\value{
+Returns data with rownames containing sequence of dates as character.
+}
+\description{
+Changes rownames to a character string containing dates. Common format of
+dates or other index is needed if the length of the \emph{instrument} and
+reduced form residuals differ. In such case, the identification is based on
+a common subset of residuals and the instrument.
+}
diff --git a/vignettes/instrument_example.Rmd b/vignettes/instrument_example.Rmd
new file mode 100644
index 0000000..9e639e1
--- /dev/null
+++ b/vignettes/instrument_example.Rmd
@@ -0,0 +1,83 @@
+---
+title: "Proxy SVAR Identification"
+output: html_document
+---
+
+```{r setup, echo=FALSE, message=FALSE}
+devtools::load_all()
+```
+
+
+This demonstration shows how to perform identification using an external instrument in an SVAR setting (Stock and Watson 2012, Mertens and Ravn 2013), which we refer to as proxy SVAR or SVAR-IV. The framework implemented here allows for different lengths of residuals and the instrument. This can be useful if the data for reduced form estimation is available for longer than the instrument used for identification. To implement this in the existing framework, we use rownames and names to subset residuals and the instrument, respectively. 
+
+Consider monthly observations of industrial production, consumer price index, unemployment rate, and 1-year treasury rate from 1979 to 2019 and an instrument for monetary policy shocks from 1991 to 2017. In this case, the proxied variable is the 1-year treasury rate (GS1). We use the instrument for monetary policy shocks from Miranda-Agrippino and Ricco (2021), extended in Degasperi and Ricco (2021).
+
+```{r}
+# Access the fred_md dataset
+data <- fred_md[,c("INDPRO","CPIAUCSL", "UNRATE", "GS1")]
+
+# Transform and subset to 1979-2019
+data <- fred_transform(data, codes = c(4, 4, 1, 1))[241:732,]
+
+# Estimate a BVAR using 12 lags, default settings and very few draws
+x <- bvar(data, lags = 12, n_draw = 600L, n_burn = 100L, verbose = FALSE)
+
+# Compute IRF using an external instrument for monetary policy shocks
+data("mpi")
+irf(x) <- irf.bvar(x, bv_irf(horizon = 24L, instrument = instrument,
+                             proxyvar = "GS1"), n_thin = 1L)
+
+```
+
+The responses can be plotted by selecting the impulse variable. Note that only shocks to proxied variable are sensibly identified. 
+
+```{r, fig.align='center'}
+margins <- c(2,10,2,10)
+plot(irf(x), vars_impulse = "GS1", mar = margins)
+```
+
+Note that rownames(data) and names(instrument) need to have the same format. This could be e.g. a simple numeric index or a specific date format. In the example above we used dates in format "YYYY-MM-DD":
+
+```{r}
+range(rownames(data))
+range(names(instrument))
+```
+
+This ensures that the correct subset of residuals can be matched with the instrument.
+
+
+Special Cases
+-------
+
+As a special case we might have data and instrument of the same length.
+```{r}
+data2 <- data[145:468,]
+
+# Estimate a BVAR using 12 lags, default settings and very few draws
+x2 <- bvar(data2, lags = 12, n_draw = 600L, n_burn = 100L, verbose = FALSE)
+
+irf(x2) <- irf.bvar(x2, bv_irf(horizon = 24L, instrument = instrument,
+                             proxyvar = "GS1"), n_thin = 1L)
+```
+Notice that the residuals are shorter than the data by the amount of the *lags* used in bvar. If the objective is to use the full sample of the instrument and the full sample of residuals, we can specify data to start *lags* observations before the instrument.
+
+```{r}
+data3 <- data[133:468,]
+
+# Estimate a BVAR using 12 lags, default settings and very few draws
+x3 <- bvar(data3, lags = 12, n_draw = 600L, n_burn = 100L, verbose = FALSE)
+
+irf(x3) <- irf.bvar(x3, bv_irf(horizon = 24L, instrument = instrument,
+                             proxyvar = "GS1"), n_thin = 1L)
+```
+
+References
+-------
+Degasperi, R. and Ricco, G. (2021) Information and policy shocks in monetary surprises. Working paper, *University of Warwick*.
+
+Mertens, K., and Ravn, M. O. (2013) The dynamic effects of personal and corporate income tax changes in the United States. *American economic review*, 103(4), 1212-1247, DOI: [10.1257/aer.103.4.1212](https://doi.org/10.1257/aer.103.4.1212).
+
+Miranda-Agrippino, S., & Ricco, G. (2021) The transmission of monetary policy shocks. *American Economic Journal: Macroeconomics*, 13(3), 74-107, DOI: [10.1257/mac.20180124](https://doi.org/10.1257/mac.20180124).
+
+Stock, J. H., and Watson, M. W. (2012) Disentangling the Channels of the 2007-2009 Recession (No. w18094). *National Bureau of Economic Research*, DOI: [10.3386/w18094](https://doi.org/10.3386/w18094).
+