diff --git a/NAMESPACE b/NAMESPACE
index 4fd1b1da..4d5ccb22 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -24,6 +24,7 @@ export(set_clinical_treatment)
 export(set_demography)
 export(set_drugs)
 export(set_equilibrium)
+export(set_housing)
 export(set_mass_pev)
 export(set_mda)
 export(set_parameter_draw)
diff --git a/R/correlation.R b/R/correlation.R
index 41ea6a82..8e5fb9ef 100644
--- a/R/correlation.R
+++ b/R/correlation.R
@@ -5,7 +5,8 @@ INTS <- c(
   'smc',
   'tbv',
   'bednets',
-  'spraying'
+  'spraying',
+  'housing'
 )
 
 #' Class: Correlation parameters
diff --git a/R/parameters.R b/R/parameters.R
index 8a96989a..8b56d8dd 100644
--- a/R/parameters.R
+++ b/R/parameters.R
@@ -149,7 +149,9 @@
 #' * k0 - proportion of females bloodfed with no net; default = 0.699
 #' * spraying - boolean for if indoor spraying is enabled; default = FALSE
 #' * phi_indoors - proportion of bites taken indoors; default = 0.90
-#'
+#' * housing - boolean for if housing improvements are enabled; default = FALSE
+#' * phi_housing - fractional increase to outdoor biting due to improved housing; default = 1 simulating no increase
+#' 
 #' treatment parameters:
 #' please set treatment parameters with the convenience functions in
 #' `drug_parameters.R`
@@ -330,6 +332,8 @@ get_parameters <- function(overrides = list()) {
     # indoor spraying
     spraying = FALSE,
     phi_indoors = .90,
+    # housing improvements
+    housing = FALSE,
     # treatment
     drug_efficacy          = numeric(0),
     drug_rel_c             = numeric(0),
diff --git a/R/processes.R b/R/processes.R
index d3b11fda..055b05af 100644
--- a/R/processes.R
+++ b/R/processes.R
@@ -235,6 +235,18 @@ create_processes <- function(
     )
   }
   
+  if (parameters$housing) {
+    processes <- c(
+      processes,
+      housing_improvement(
+        variables,
+        parameters,
+        correlations
+      ),
+      house_usage_renderer(variables$house_time, renderer)
+    )
+  }
+  
   # ======================
   # Progress bar process
   # ======================
diff --git a/R/variables.R b/R/variables.R
index 0931b3f7..3481d4e9 100644
--- a/R/variables.R
+++ b/R/variables.R
@@ -26,6 +26,7 @@
 #' haven't been any
 #' * net_time - The timestep when a net was last put up (-1 if never)
 #' * spray_time - The timestep when the house was last sprayed (-1 if never)
+#' * house_time - The timestep for adaptation to the house to reduce entry/kill vectors (-1 if never)
 #' * infectivity - The onward infectiousness to mosquitos
 #' * drug - The last prescribed drug
 #' * drug_time - The timestep of the last drug
@@ -196,6 +197,7 @@ create_variables <- function(parameters) {
   # Init vector controls
   net_time <- individual::IntegerVariable$new(rep(-1, size))
   spray_time <- individual::IntegerVariable$new(rep(-1, size))
+  house_time <- individual::IntegerVariable$new(rep(-1, size))
 
   variables <- list(
     state = state,
@@ -219,7 +221,8 @@ create_variables <- function(parameters) {
     pev_profile = pev_profile,
     tbv_vaccinated = tbv_vaccinated,
     net_time = net_time,
-    spray_time = spray_time
+    spray_time = spray_time,
+    house_time = house_time
   )
 
   # Add variables for individual mosquitoes
diff --git a/R/vector_control.R b/R/vector_control.R
index 29c92d35..15946f47 100644
--- a/R/vector_control.R
+++ b/R/vector_control.R
@@ -11,7 +11,7 @@ prob_bitten <- function(
   parameters
   ) {
   n <- parameters$human_population
-  if (!(parameters$bednets || parameters$spraying)) {
+  if (!(parameters$bednets || parameters$spraying || parameters$housing)) {
     return(
       list(
         prob_bitten_survives = rep(1, n),
@@ -65,6 +65,7 @@ prob_bitten <- function(
       js_prime,
       parameters$k0
     )
+    spray_on = 1
     rs_comp <- 1 - rs
     ss <- rep(1, n)
     ss[protected_index] <- prob_survives_spraying(
@@ -72,30 +73,76 @@ prob_bitten <- function(
       parameters$k0
     )
   } else {
-    phi_indoors <- 0
+    spray_on = 0
     rs <- 0
     rs_comp <- 1
     ss <- 1
   }
+  
+  if (parameters$housing) {
+    phi_housing <- parameters$phi_housing[[species]]  ## if a change in behaviour is caused by the housing changes, we can increase outdoor biting proportion
+    phi_indoors <- parameters$phi_indoors[[species]]
+    house_time <- variables$house_time$get_values()
+    since_housing <- timestep - house_time
+    matches <- match(house_time, parameters$housing_timesteps)
+    rh <- prob_repelled_house(matches, since_housing, species, parameters) ## if housing prevents entry to house, we increase proportion needing to repeat foraging
+    sh <- prob_survives_house(rh, matches, since_housing, species, parameters)
+    unused <- house_time == -1
+    sh[unused] <- 1
+    rh[unused] <- 0
+  } else {
+    phi_housing <- 1
+    rh <- 0
+    sh <- 1
+  }
 
+  if ((!parameters$housing & !parameters$spraying)) {
+    phi_indoors <- 0 ## we want phi_indoors to be applied if housing is on
+  }
+  
   list(
     prob_bitten_survives = (
-      1 - phi_indoors +
-      phi_bednets * rs_comp * sn * ss +
-      (phi_indoors - phi_bednets) * rs_comp * ss
+      1 - phi_indoors * phi_housing + ##
+        (phi_bednets * phi_housing * rs_comp * sn * ss * sh^2) +  ## * sh if some mortality from housing
+        ((phi_indoors * phi_housing - phi_bednets * phi_housing) * rs_comp * ss * sh^2) ## * sh if some mortality from housing
     ),
     prob_bitten = (
-      1 - phi_indoors +
-      phi_bednets * rs_comp * sn +
-      (phi_indoors - phi_bednets) * rs_comp
+      1 - phi_indoors * phi_housing +
+        (phi_bednets * phi_housing * rs_comp * sn * sh) + ## * sh if some mortality from housing
+        ((phi_indoors * phi_housing - phi_bednets * phi_housing) * rs_comp * sh) ## * sh if some mortality from housing
     ),
     prob_repelled = (
-      phi_bednets * rs_comp * rn +
-      phi_indoors * rs
+      phi_bednets * phi_housing * rs_comp * sh * rn * sh + 
+        phi_indoors * phi_housing * rs * sh * sh * spray_on + 
+        phi_indoors * phi_housing * rh
     )
   )
 }
 
+#' @title Simulate housing improvements
+#' @description simulates improved housing so that harder for vectors to get indoors
+#' from `set_housing` and correlation parameters from
+#' `get_correlation_parameters`
+#'
+#' @param variables list of variables in the model
+#' @param parameters the model parameters
+#' @param correlations correlation parameters
+#' @noRd
+housing_improvement <- function(variables, parameters, correlations) {
+  function(timestep) {
+    matches <- timestep == parameters$housing_timesteps
+    if (any(matches)) {
+      target <- which(sample_intervention(
+        seq(parameters$human_population),
+        'housing',
+        parameters$housing_coverages[matches],
+        correlations
+      ))
+      variables$house_time$queue_update(timestep, target)
+    }
+  }
+}
+
 #' @title Indoor spraying
 #' @description models indoor residual spraying according to the strategy
 #' from `set_spraying` and correlation parameters from
@@ -106,7 +153,7 @@ prob_bitten <- function(
 #' @param correlations correlation parameters
 #' @noRd
 indoor_spraying <- function(spray_time, parameters, correlations) {
-  function(timestep) {
+      function(timestep) {
     matches <- timestep == parameters$spraying_timesteps
     if (any(matches)) {
       target <- which(sample_intervention(
@@ -187,6 +234,22 @@ spraying_decay <- function(t, theta, gamma) {
   1 / (1 + exp(-(theta + gamma * t)))
 }
 
+house_decay <- function(t, gamma) {
+  exp(-t / gamma)
+}
+
+prob_repelled_house <- function(matches, dt, species, parameters) {
+  rhm <- parameters$house_rhm[matches, species]
+  gammah <- parameters$house_gammah[matches]
+  (parameters$house_rh[matches, species] - rhm) * house_decay(dt, gammah) + rhm
+}
+
+prob_survives_house <- function(rh, matches, dt, species, parameters) {
+  dh0 <- parameters$house_dh0[matches, species]
+  dh <- dh0 * house_decay(dt, parameters$house_gammah[matches])
+  1 - rh - dh
+}
+
 net_usage_renderer <- function(net_time, renderer) {
   function(t) {
     renderer$render(
@@ -196,3 +259,13 @@ net_usage_renderer <- function(net_time, renderer) {
     )
   }
 }
+
+house_usage_renderer <- function(house_time, renderer) {
+  function(t) {
+    renderer$render(
+      'n_use_house_adaptation',
+      house_time$get_index_of(-1)$not(TRUE)$size(),
+      t
+    )
+  }
+}
diff --git a/R/vector_control_parameters.R b/R/vector_control_parameters.R
index 97c3ce56..2602c131 100644
--- a/R/vector_control_parameters.R
+++ b/R/vector_control_parameters.R
@@ -179,3 +179,61 @@ get_init_carrying_capacity <- function(parameters){
   names(init_cc) <- parameters$species
   return(init_cc)
 }
+
+#' @title Parameterise a housing improvement strategy
+#'
+#' @description The model will simulate improved housing at `timesteps` to a random
+#' sample of the entire human population. The sample size will be a proportion
+#' of the human population taken from the corresponding `coverages`.
+#' The sample _can_ contain humans who have already benefited from housing.
+#'
+#' If a human in the sample lives in a house that has been improved to reduce
+#' mosquito entry, the efficacy of the
+#' housing improvement will be high - as determined by the parameter rn_house
+#'
+#' The structure for the housing model will be documented in a publication
+#' Sherrard-Smith et al in prep
+#'
+#' @param parameters a list of parameters to modify
+#' @param timesteps the timesteps at which to distribute housing adaptations
+#' @param coverages the proportion of the human population who reside in protected housing
+#' @param dh0 a matrix of death probabilities for each species over time.
+#' With nrows=length(timesteps), ncols=length(species)
+#' @param rh a matrix of repelling probabilities for each species over time
+#' With nrows=length(timesteps), ncols=length(species)
+#' @param rhm a matrix of minimum repelling probabilities for each species over time
+#' With nrows=length(timesteps), ncols=length(species)
+#' @param gammah a vector of house adaptation half-lives for each distribution timestep
+#' @export
+set_housing <- function(
+    parameters,
+    timesteps,
+    coverages,
+    phi_housing,
+    dh0,
+    rh,
+    rhm,
+    gammah
+) {
+  stopifnot(all(coverages >= 0) && all(coverages <= 1))
+  lengths <- vnapply(list(coverages, gammah), length)
+  if (!all(lengths == length(timesteps))) {
+    stop('timesteps and time-varying parameters must align')
+  }
+  for (x in list(dh0, rh, rhm)) {
+    if (ncol(x) != length(parameters$species)) {
+      stop('death and repelling probabilities rows need to align with species')
+    }
+    if (nrow(x) != length(timesteps)) {
+      stop('death and repelling probabilities columns need to align with timesteps')
+    }
+  }
+  parameters$housing <- TRUE
+  parameters$housing_timesteps <- timesteps
+  parameters$housing_coverages <- coverages
+  parameters$house_dh0 <- dh0
+  parameters$house_rh <- rh
+  parameters$house_rhm <- rhm
+  parameters$house_gammah <- gammah
+  parameters
+}
\ No newline at end of file
diff --git a/R/vector_parameters.R b/R/vector_parameters.R
index 2024d253..267b25ac 100644
--- a/R/vector_parameters.R
+++ b/R/vector_parameters.R
@@ -6,6 +6,7 @@
 #' Q0: 0.92
 #' phi_bednets: 0.85
 #' phi_indoors: 0.90
+#' phi_housing: 1
 #' mum: 0.132
 #'
 #' parameters from:
@@ -18,6 +19,7 @@ gamb_params <- list(
   Q0 = .92,
   phi_bednets = .85,
   phi_indoors = .90,
+  phi_housing = 1,
   mum = .132
 )
 
@@ -29,6 +31,7 @@ gamb_params <- list(
 #' Q0: 0.71
 #' phi_bednets: 0.8
 #' phi_indoors: 0.86
+#' phi_housing: 1
 #' mum: 0.132
 #'
 #' parameters from:
@@ -41,6 +44,7 @@ arab_params <- list(
   Q0 = .71,
   phi_bednets = .8,
   phi_indoors = .86,
+  phi_housing = 1,
   mum = .132
 )
 
@@ -52,6 +56,7 @@ arab_params <- list(
 #' Q0: 0.94
 #' phi_bednets: 0.78
 #' phi_indoors: 0.87
+#' phi_housing: 1
 #' mum: 0.112
 #'
 #' parameters from:
@@ -64,6 +69,7 @@ fun_params <- list(
   Q0 = .94,
   phi_bednets = .78,
   phi_indoors = .87,
+  phi_housing = 1,
   mum = .112
 )
 
@@ -75,6 +81,7 @@ fun_params <- list(
 #' Q0: 0.21
 #' phi_bednets: 0.57
 #' phi_indoors: 0.37
+#' phi_housing: 1
 #' mum: 0.112
 #'
 #' parameters reference:
@@ -91,6 +98,7 @@ steph_params <- list(
   Q0 = 0.21,
   phi_bednets = 0.52186,
   phi_indoors = 0.4776,
+  phi_housing = 1,
   mum = .112
 )
 
@@ -114,6 +122,7 @@ set_species <- function(parameters, species, proportions) {
     'Q0',
     'phi_bednets',
     'phi_indoors',
+    'phi_housing',
     'mum'
   )
   for (key in keys) {
diff --git a/man/arab_params.Rd b/man/arab_params.Rd
index 144bbf70..c8d71929 100644
--- a/man/arab_params.Rd
+++ b/man/arab_params.Rd
@@ -5,7 +5,7 @@
 \alias{arab_params}
 \title{Preset parameters for the An. arabiensis vector}
 \format{
-An object of class \code{list} of length 7.
+An object of class \code{list} of length 8.
 }
 \usage{
 arab_params
@@ -21,6 +21,7 @@ foraging_time: .69
 Q0: 0.71
 phi_bednets: 0.8
 phi_indoors: 0.86
+phi_housing: 1
 mum: 0.132
 
 parameters from:
diff --git a/man/fun_params.Rd b/man/fun_params.Rd
index 6063bfbc..a53e1f29 100644
--- a/man/fun_params.Rd
+++ b/man/fun_params.Rd
@@ -5,7 +5,7 @@
 \alias{fun_params}
 \title{Preset parameters for the An. funestus vector}
 \format{
-An object of class \code{list} of length 7.
+An object of class \code{list} of length 8.
 }
 \usage{
 fun_params
@@ -21,6 +21,7 @@ foraging_time: .69
 Q0: 0.94
 phi_bednets: 0.78
 phi_indoors: 0.87
+phi_housing: 1
 mum: 0.112
 
 parameters from:
diff --git a/man/gamb_params.Rd b/man/gamb_params.Rd
index fef4b29f..0af77034 100644
--- a/man/gamb_params.Rd
+++ b/man/gamb_params.Rd
@@ -5,7 +5,7 @@
 \alias{gamb_params}
 \title{Preset parameters for the An. gambiae s.s vector}
 \format{
-An object of class \code{list} of length 7.
+An object of class \code{list} of length 8.
 }
 \usage{
 gamb_params
@@ -21,6 +21,7 @@ foraging_time: .69
 Q0: 0.92
 phi_bednets: 0.85
 phi_indoors: 0.90
+phi_housing: 1
 mum: 0.132
 
 parameters from:
diff --git a/man/get_parameters.Rd b/man/get_parameters.Rd
index 54188c01..b9b01ad7 100644
--- a/man/get_parameters.Rd
+++ b/man/get_parameters.Rd
@@ -166,6 +166,8 @@ please set vector control strategies using \code{set_betnets} and \code{set_spra
 \item k0 - proportion of females bloodfed with no net; default = 0.699
 \item spraying - boolean for if indoor spraying is enabled; default = FALSE
 \item phi_indoors - proportion of bites taken indoors; default = 0.90
+\item housing - boolean for if housing improvements are enabled; default = FALSE
+\item phi_housing - fractional increase to outdoor biting due to improved housing; default = 1 simulating no increase
 }
 
 treatment parameters:
diff --git a/man/set_housing.Rd b/man/set_housing.Rd
new file mode 100644
index 00000000..c37d68a9
--- /dev/null
+++ b/man/set_housing.Rd
@@ -0,0 +1,48 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/vector_control_parameters.R
+\name{set_housing}
+\alias{set_housing}
+\title{Parameterise a housing improvement strategy}
+\usage{
+set_housing(
+  parameters,
+  timesteps,
+  coverages,
+  phi_housing,
+  dh0,
+  rh,
+  rhm,
+  gammah
+)
+}
+\arguments{
+\item{parameters}{a list of parameters to modify}
+
+\item{timesteps}{the timesteps at which to distribute housing adaptations}
+
+\item{coverages}{the proportion of the human population who reside in protected housing}
+
+\item{dh0}{a matrix of death probabilities for each species over time.
+With nrows=length(timesteps), ncols=length(species)}
+
+\item{rh}{a matrix of repelling probabilities for each species over time
+With nrows=length(timesteps), ncols=length(species)}
+
+\item{rhm}{a matrix of minimum repelling probabilities for each species over time
+With nrows=length(timesteps), ncols=length(species)}
+
+\item{gammah}{a vector of house adaptation half-lives for each distribution timestep}
+}
+\description{
+The model will simulate improved housing at \code{timesteps} to a random
+sample of the entire human population. The sample size will be a proportion
+of the human population taken from the corresponding \code{coverages}.
+The sample \emph{can} contain humans who have already benefited from housing.
+
+If a human in the sample lives in a house that has been improved to reduce
+mosquito entry, the efficacy of the
+housing improvement will be high - as determined by the parameter rn_house
+
+The structure for the housing model will be documented in a publication
+Sherrard-Smith et al in prep
+}
diff --git a/man/steph_params.Rd b/man/steph_params.Rd
index 99b81bcc..129e092e 100644
--- a/man/steph_params.Rd
+++ b/man/steph_params.Rd
@@ -5,7 +5,7 @@
 \alias{steph_params}
 \title{Preset parameters for the An. stephensi vector}
 \format{
-An object of class \code{list} of length 7.
+An object of class \code{list} of length 8.
 }
 \usage{
 steph_params
@@ -21,6 +21,7 @@ foraging_time: .69
 Q0: 0.21
 phi_bednets: 0.57
 phi_indoors: 0.37
+phi_housing: 1
 mum: 0.112
 
 parameters reference:
diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp
index f5c226fd..affb233d 100644
--- a/src/RcppExports.cpp
+++ b/src/RcppExports.cpp
@@ -260,7 +260,7 @@ BEGIN_RCPP
 END_RCPP
 }
 
-RcppExport SEXP run_testthat_tests();
+RcppExport SEXP run_testthat_tests(void);
 
 static const R_CallMethodDef CallEntries[] = {
     {"_malariasimulation_create_adult_mosquito_model", (DL_FUNC) &_malariasimulation_create_adult_mosquito_model, 5},
diff --git a/vignettes/VectorControl_Housing.Rmd b/vignettes/VectorControl_Housing.Rmd
new file mode 100644
index 00000000..a88d9354
--- /dev/null
+++ b/vignettes/VectorControl_Housing.Rmd
@@ -0,0 +1,125 @@
+---
+title: "Vector Control: Housing adaptations"
+output: 
+  rmarkdown::html_vignette:
+vignette: >
+  %\VignetteIndexEntry{Vector Control: Housing improvements}
+  %\VignetteEngine{knitr::rmarkdown}
+  %\VignetteEncoding{UTF-8}
+---
+
+```{r, include = FALSE}
+knitr::opts_chunk$set(
+  collapse = TRUE,
+  comment = "#>",
+  dpi=300,
+  fig.width=7
+)
+```
+
+```{r setup}
+# Load the requisite packages:
+library(malariasimulation)
+# Set colour palette:
+cols <- c("#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7")
+```
+
+Housing adaptations are continually occurring in malaria-endemic parts of the world. Improvements could reduce the penetrability of housing to vectors. This will be a longer-lasting protective barrier effect than other 'indoor interventions' like ITNs or spraying. The model provides the flexibility to turn on 'improved' housing using a mechanism that adjusts the mosquito foraging pathway. The female mosquito now has some probability of being repelled by the adjustments to the house prior to any effects of spray or ITNs. The impacts from the three types of interventions are coded to be multiplicative and independent. When housing is 0, there is no impact. The model provides the user with the flexibility to specify parameters that describe the housing campaign (e.g. timing, coverage, and efficacy). It is currently assumed that housing improvements will represent a permanent change. We can illustrate this below:
+
+We can create a few plotting functions to visualise the output.
+```{r}
+# Plotting functions
+plot_prev <- function() {
+  plot(x = output$timestep, y = output$n_detect_730_3650 / output$n_730_3650, 
+       type = "l", col = cols[3], lwd = 1,
+       xlab = "Time (days)", ylab = expression(paste(italic(Pf),"PR"[2-10])),
+       xaxs = "i", yaxs = "i", ylim = c(0, 1))
+  lines(x = output_control$timestep, y = output_control$n_detect_730_3650 / output_control$n_730_3650,
+        col = cols[5], lwd = 1)
+  lines(x = output_nets_housing$timestep, y = output_nets_housing$n_detect_730_3650 / output_nets_housing$n_730_3650,
+        col = cols[1], lwd = 1)
+  abline(v = bednetstimesteps, col = "black", lty = 2, lwd = 1)
+  abline(v = housingtimesteps, col = "darkred", lty = 2, lwd = 1)
+  
+  text(x = bednetstimesteps + 10, y = 0.95, labels = "Bed net int.", adj = 0, cex = 0.8)
+  text(x = housingtimesteps + 10, y = 0.95, labels = "Housing int.", adj = 0, cex = 0.8)
+  grid(lty = 2, col = "grey80", lwd = 0.5)
+  legend("bottomleft", box.lty = 0, bg = "white",
+         legend = c("Prevalence for housing scenario",
+                    "Prevalence for nets and housing scenario",
+                    "Prevalence for control scenario"),
+         col = c(cols[3], cols[4],cols[1],cols[5]), lty = c(1,1,1,1), lwd = 2, cex = 0.8, y.intersp = 1.3)
+}
+```
+  
+## Setting housing parameters
+
+### Parameterisation
+
+Use the `get_parameters()` function to generate the list of parameters for a perennial profile, accepting the default values to run the simulation from an equilibrium starting point. 
+
+```{r}
+year <- 365
+sim_length <- 6 * year
+human_population <- 1000
+starting_EIR <- 50
+
+simparams <- get_parameters(
+  list(human_population = human_population)
+)
+
+simparams <- set_equilibrium(parameters = simparams, init_EIR = starting_EIR)
+
+output_control <- run_simulation(timesteps = sim_length, parameters = simparams)
+```
+
+#### A note on mosquito species  
+It is also possible to use the `set_species()` function to account for different mosquito species in the simulation. In this case, the matrices would need to have additional columns corresponding to each mosquito species. If you are not already familiar with the `set_species()` function, see the [Mosquito Species](https://mrc-ide.github.io/malariasimulation/articles/SetSpecies.html) vignette.  
+
+The default parameters are set to model *Anopheles gambiae*. 
+```{r}
+simparams$species
+simparams$species_proportions
+```
+  
+### Simulation
+
+Having established a base set of parameters, we can now create a copy of this parameter list and update it to specify that housing has been altered in some way to reduce access to mosquitoes. In the example below, we simulate that 30% of housing is 'improved' in the 2nd year, and that nets are used by a random 50% of the population every three years. We can change the 'repellence' introduced from households multiple times and differently for each mosquito species potentially.
+  
+```{r, fig.align = 'center', out.width='100%'}
+housingparams <- set_housing(
+  simparams,
+  timesteps = housingtimesteps,
+  coverages = c(0.8,0.8,0.9),
+  phi_housing = c(1,1,1),
+  dh0 = matrix(c(.2, .2), nrow = 3, ncol = 2), 
+  rh = matrix(c(.5, .5), nrow = 3, ncol = 2), 
+  rhm = matrix(c(.3, .3), nrow = 3, ncol = 2), 
+  gammah = rep(0.5 * 365, 3) 
+)
+output <- run_simulation(timesteps = sim_length, parameters = housingparams)
+
+
+bednetstimesteps <- c(1, 4) * year # The bed nets will be distributed at the end of the first and the 4th year. 
+
+bednetparams <- set_bednets(
+  housingparams,
+  timesteps = bednetstimesteps,
+  coverages = c(.5, .8),  
+  retention = 5 * year, 
+  dn0 = matrix(c(.533, .533), nrow = 2, ncol = 1), 
+  rn = matrix(c(.56, .56), nrow = 2, ncol = 1), 
+  rnm = matrix(c(.24, .24), nrow = 2, ncol = 1), 
+  gamman = rep(2.64 * 365, 2) 
+)
+
+output_nets_housing <- run_simulation(timesteps = sim_length, parameters = bednetparams)
+```
+
+### Visualisation
+
+```{r, fig.align = 'center', out.width='100%'}
+# Plot prevalence
+plot_prev()
+```
+