The goal of this assignment is evaluate the Watts-Strogatz model (WS) by running Watts & Strogatzâ€™ (1998) key simulation showing the impact of p on clustering and shortest path lengths.

This assignment contains of 4 steps.

- Place your model inside a function.
- Running the simulation.
- Controlling for the results of a regular lattice.
- Plotting the results.

- The first step of running the simulation is to put the Watts-Strogatz model into its own function, which is actually easier than it sounds. Functions are created in R using the code below.

```
function_name <- function(ARGUMENT1, ARGUMENT2, ARGUMENT3){
CODE_HERE
}
```

- What you need to do is three things. First, place your Watts & Strogatz model in-between the curly brackets (inplace of
`CODE_HERE`

). If you donâ€™t have a running version of the model, use this one:

```
# get edges
edges = matrix(ncol = 3, nrow = size * max_dist)
ind = 0
for(i in 1:(size-1)){
for(j in (i+1):size){
dist = min(c(abs(j - i),abs((j - size)-i)))
if(i != j & dist <= max_dist){
ind = ind + 1
edges[ind, ] = c(i, j, dist)
}
}
}
# sort edges by distance
edges = edges[order(edges[,3],edges[,1],edges[,2]),]
# fill matrix
network = matrix(0, ncol = size, nrow = size)
network[edges[,1:2]] = network[edges[,2:1]] = 1
# rewire
for(i in 1:nrow(edges)){
if(runif(1,0,1) < p){
new_end = sample(1:size,1)
if(new_end != edges[i, 1] & network[edges[i,1], new_end] == 0){
network[edges[i,1], edges[i,2]] = network[edges[i,2], edges[i,1]] = 0
network[edges[i,1], new_end] = network[new_end, edges[i,1]] = 1
}
}
}
network
}
```

Now, change the names of

`ARGUMENT1`

,`ARGUMENT2`

, and`ARGUMENT3`

to the names of the three objects that control the size of the network, the neighborhood size (i.e.,`max_dist`

), and the probability of rewiring (in that order). Finally, change the name of the function to`watts_strogatz`

. Thatâ€™s it. Run the code and you can use your function like any other.Try running

`watts_strogatz(10, 2, .2)`

. Does it work? Does the resulting matrix appear correct? If not, try to find the error before you proceed.

In their simulation, Watts and Strogatz varied the probability of rewiring edges to evaluate its impact on clustering and average shortest path length. To do the same, we need to find a way to run the same model for different values of the rewiring probability. Solution: a for-loop. We begin by defining values for the rewiring probability using

`ps <- seq(start, end, step)`

where`start`

is the smallest value,`end`

the largest, and`step`

the step size of the sequence. Choose`start`

and`end`

to be`.00001`

and`1`

, and step to be some number that produces, for now, no more than say 10 or 20 values.Now create a loop that iterates over those

`ps`

and another one inside that loop that repeats the simulation for each value in`ps`

exactly`n`

times. Just like the code below. Chose`n`

to be a small number, say 10, for now. Later you can increase it to increase the precision of your simulation. Alright, this is the backbone of our simulation.

```
n = 10
for(i in 1:length(ps)){
for(j in 1:n){
CODE_HERE
}
}
```

Letâ€™s figure out how to calculate clustering and shortest path lengths. To do this, we want to use the

`igraph`

package. Once you installed the package (`install.packages('igraph')`

) and loaded it (`library(igraph)`

), you can use`graph_from_adjacency_matrix()`

to load the matrix produced by your`watts_strogatz`

-function into the`igraph`

format. That is, assume`network`

is an object created from`watts_strogatz`

, use`graph_from_adjacency_matrix(network, mode = 'undirected')`

and assign the result to a new object called`igraph_network`

.Now you can use the igraph functions

`transitivity()`

and`average.path.length()`

on the new`igraph_network`

-object to compute clustering and average shortest path length. In using`transitivity`

, make sure to set`type`

to`"localaverage"`

to compute the type of clustering that Watts & Strogatz used.Ok, now you pretty much got all the ingredients. Whatâ€™s left is storing. Create a 2-column matrix with

`n * length(ps)`

rows called`result`

. Also create an`index`

object and set it to`1`

. Create both before the loops. Then inside the inner loop assign the result of`transitivity()`

to`result[index, 1]`

and the result of`average.path.length()`

to`result[index, 2]`

.Now, you are technically ready to run the simulation. Before starting it, however, I recommend placing an

`print(i)`

between the inner and outer loop. This will help you keep track of the progress of the loop. Now run the simulation once with network`size`

equal 10 and`max_dist`

equal 2. If it takes to long, reduce the length of`ps`

and/or the value of`n`

. If it runs fine, set the values of`size`

and`max_dist`

to what Watts & Strogatz used in their simualtion.

If you pay close attention, you will see that Watts and Strogatz do not show the raw clustering coefficients and average shortest path length, but those normalized for the results of regular lattices. Obtain the necessary normalization constants, by creating networks using

`watts_strogatz(size, max_dist, p = 0)`

- remember`p = 0`

means regular lattice - and evaluating its clustering and shortest path length using the same functions as before. Store the results as`C0`

and`L0`

, respectively.Now normalize the simulation results using

`result[, 1] <- result[, 1] / C0`

and`result[, 2] <- result[, 2] / L0`

.

- Ok, you ran the simulation, got results, normalized them, and now its finally time to inspect the results. First, aggregate the results of the repeated simulations. Do this using
`agg_C <- tapply(res[,1], rep(ps, rep(n, length(ps))), mean)`

and, respectively, the same for the average shortest path lengths. Now, plot the relationship between`ps`

and`agg_C`

using`plot(ps, agg_C, type = 'l', ylim = c(0, 1), log = 'x')`

. Then run`lines(ps, agg_L)`

, which will add a dashed line for the average shortest path lengths. What do you see? Does it look the same as in Watts & Strogatz (1998)?

```
require(igraph)
watts_strogatz = function(size = 10, max_dist = 2, p = 0, igraph = TRUE){
# get edges
edges = matrix(ncol = 3, nrow = size * max_dist)
ind = 0
for(i in 1:(size-1)){
for(j in (i+1):size){
dist = min(c(abs(j - i),abs((j - size)-i)))
if(i != j & dist <= max_dist){
ind = ind + 1
edges[ind, ] = c(i, j, dist)
}
}
}
# sort edges by distance
edges = edges[order(edges[,3],edges[,1],edges[,2]),]
# fill matrix
network = matrix(0, ncol = size, nrow = size)
network[edges[,1:2]] = network[edges[,2:1]] = 1
# rewire
for(i in 1:nrow(edges)){
if(runif(1,0,1) < p){
new_end = sample(1:size,1)
if(new_end != edges[i, 1] & network[edges[i,1], new_end] == 0){
network[edges[i,1], edges[i,2]] = network[edges[i,2], edges[i,1]] = 0
network[edges[i,1], new_end] = network[new_end, edges[i,1]] = 1
}
}
}
network
}
# simulate
ps <- seq(.00001, 1, length = 20)
n = 10
res = matrix(ncol = 2, nrow = length(ps) * n)
ind = 0
for(i in 1:length(ps)){
print(i)
for(j in 1:n){
ind = ind + 1
network <- watts_strogatz(1000, 10, ps[i])
igraph_network = graph_from_adjacency_matrix(network, mode = "undirected")
res[ind, 1] = transitivity(igraph_network, type = 'localaverage')
res[ind, 2] = average.path.length(igraph_network)
}
}
```

```
## [1] 1
## [1] 2
## [1] 3
## [1] 4
## [1] 5
## [1] 6
## [1] 7
## [1] 8
## [1] 9
## [1] 10
## [1] 11
## [1] 12
## [1] 13
## [1] 14
## [1] 15
## [1] 16
## [1] 17
## [1] 18
## [1] 19
## [1] 20
```

```
# normalize
network_0 = graph_from_adjacency_matrix(watts_strogatz(1000, 10, 0), mode = "undirected")
res[,1] = res[,1] / transitivity(network_0, type = 'localaverage')
res[,2] = res[,2] / average.path.length(network_0)
# aggregate & plot
agg_C = tapply(res[,1], rep(ps, rep(n, length(ps))), mean)
agg_L = tapply(res[,2], rep(ps, rep(n, length(ps))), mean)
# plot
plot(ps, agg_C, type = 'l', ylim = c(0, 1), log = 'x')
lines(ps, agg_L, lty = 2)
```