dat2=NULL

#Gore Bay
dat1=read.csv("mm6092920.txt",header=FALSE,skip=4,,col.names=c("Year","Jan","Jan Status","Feb","Feb Status","Mar","Mar Status","Apr","Apr Status","May","May Status","Jun","Jun Status","Jul","Jul Status","Aug","Aug Status","Sept","Sept Status","Oct","Oct Status","Nov","Nov Status","Dec","Dec Status","Annual","Annual Status","Spring","Spring Status","Summer","Summer Status","Fall","Fall Status","Winter","Winter Status"))

dat2=cbind(dat1$Year,dat1$Jan,dat1$Feb,dat1$Mar,dat1$Apr,dat1$May,dat1$Jun,dat1$Jul,dat1$Aug,dat1$Sept,dat1$Oct,dat1$Nov,dat1$Dec,rep(45.88,length(dat1$Year)),rep(-82.57,length(dat1$Year)),rep(194,length(dat1$Year)))
colnames(dat2)=c("Year","Jan","Feb","Mar","Apr","May","Jun","Jul","Aug","Sept","Oct","Nov","Dec","latitude","longitude","elevation")

#BIG TROUT LAKE
#dat1=read.csv("mm6010735.txt",header=FALSE,skip=4,,col.names=c("Year","Jan","Jan Status","Feb","Feb Status","Mar","Mar Status","Apr","Apr Status","May","May Status","Jun","Jun Status","Jul","Jul Status","Aug","Aug Status","Sept","Sept Status","Oct","Oct Status","Nov","Nov Status","Dec","Dec Status","Annual","Annual Status","Spring","Spring Status","Summer","Summer Status","Fall","Fall Status","Winter","Winter Status"))

#dat2=rbind(dat2,cbind(dat1$Year,dat1$Jan,dat1$Feb,dat1$Mar,dat1$Apr,dat1$May,dat1$Jun,dat1$Jul,dat1$Aug,dat1$Sept,dat1$Oct,dat1$Nov,dat1$Dec,rep(53.83,length(dat1$Year)),rep(-89.87,length(dat1$Year)),rep(224,length(dat1$Year))))

#LANSDOWNE HOUSE
dat1=read.csv("mm6014353.txt",header=FALSE,skip=4,,col.names=c("Year","Jan","Jan Status","Feb","Feb Status","Mar","Mar Status","Apr","Apr Status","May","May Status","Jun","Jun Status","Jul","Jul Status","Aug","Aug Status","Sept","Sept Status","Oct","Oct Status","Nov","Nov Status","Dec","Dec Status","Annual","Annual Status","Spring","Spring Status","Summer","Summer Status","Fall","Fall Status","Winter","Winter Status"))

dat2=rbind(dat2,cbind(dat1$Year,dat1$Jan,dat1$Feb,dat1$Mar,dat1$Apr,dat1$May,dat1$Jun,dat1$Jul,dat1$Aug,dat1$Sept,dat1$Oct,dat1$Nov,dat1$Dec,rep(52.23,length(dat1$Year)),rep(-87.88,length(dat1$Year)),rep(255,length(dat1$Year))))

#PICKLE LAKE
dat1=read.csv("mm6016529.txt",header=FALSE,skip=4,,col.names=c("Year","Jan","Jan Status","Feb","Feb Status","Mar","Mar Status","Apr","Apr Status","May","May Status","Jun","Jun Status","Jul","Jul Status","Aug","Aug Status","Sept","Sept Status","Oct","Oct Status","Nov","Nov Status","Dec","Dec Status","Annual","Annual Status","Spring","Spring Status","Summer","Summer Status","Fall","Fall Status","Winter","Winter Status"))

dat2=rbind(dat2,cbind(dat1$Year,dat1$Jan,dat1$Feb,dat1$Mar,dat1$Apr,dat1$May,dat1$Jun,dat1$Jul,dat1$Aug,dat1$Sept,dat1$Oct,dat1$Nov,dat1$Dec,rep(51.45,length(dat1$Year)),rep(-90.22,length(dat1$Year)),rep(386,length(dat1$Year))))


#RED LAKE
dat1=read.csv("mm6016970.txt",header=FALSE,skip=4,,col.names=c("Year","Jan","Jan Status","Feb","Feb Status","Mar","Mar Status","Apr","Apr Status","May","May Status","Jun","Jun Status","Jul","Jul Status","Aug","Aug Status","Sept","Sept Status","Oct","Oct Status","Nov","Nov Status","Dec","Dec Status","Annual","Annual Status","Spring","Spring Status","Summer","Summer Status","Fall","Fall Status","Winter","Winter Status"))

dat2=rbind(dat2,cbind(dat1$Year,dat1$Jan,dat1$Feb,dat1$Mar,dat1$Apr,dat1$May,dat1$Jun,dat1$Jul,dat1$Aug,dat1$Sept,dat1$Oct,dat1$Nov,dat1$Dec,rep(51.07,length(dat1$Year)),rep(-93.8,length(dat1$Year)),rep(386,length(dat1$Year))))

#FORT FRANCES
dat1=read.csv("mm6022474.txt",header=FALSE,skip=4,,col.names=c("Year","Jan","Jan Status","Feb","Feb Status","Mar","Mar Status","Apr","Apr Status","May","May Status","Jun","Jun Status","Jul","Jul Status","Aug","Aug Status","Sept","Sept Status","Oct","Oct Status","Nov","Nov Status","Dec","Dec Status","Annual","Annual Status","Spring","Spring Status","Summer","Summer Status","Fall","Fall Status","Winter","Winter Status"))

dat2=rbind(dat2,cbind(dat1$Year,dat1$Jan,dat1$Feb,dat1$Mar,dat1$Apr,dat1$May,dat1$Jun,dat1$Jul,dat1$Aug,dat1$Sept,dat1$Oct,dat1$Nov,dat1$Dec,rep(48.65,length(dat1$Year)),rep(-93.43,length(dat1$Year)),rep(342,length(dat1$Year))))

#MINE CENTRE
dat1=read.csv("mm6025205.txt",header=FALSE,skip=4,,col.names=c("Year","Jan","Jan Status","Feb","Feb Status","Mar","Mar Status","Apr","Apr Status","May","May Status","Jun","Jun Status","Jul","Jul Status","Aug","Aug Status","Sept","Sept Status","Oct","Oct Status","Nov","Nov Status","Dec","Dec Status","Annual","Annual Status","Spring","Spring Status","Summer","Summer Status","Fall","Fall Status","Winter","Winter Status"))

dat2=rbind(dat2,cbind(dat1$Year,dat1$Jan,dat1$Feb,dat1$Mar,dat1$Apr,dat1$May,dat1$Jun,dat1$Jul,dat1$Aug,dat1$Sept,dat1$Oct,dat1$Nov,dat1$Dec,rep(48.8,length(dat1$Year)),rep(-92.6,length(dat1$Year)),rep(361,length(dat1$Year))))

#DRYDEN
dat1=read.csv("mm6032125.txt",header=FALSE,skip=4,,col.names=c("Year","Jan","Jan Status","Feb","Feb Status","Mar","Mar Status","Apr","Apr Status","May","May Status","Jun","Jun Status","Jul","Jul Status","Aug","Aug Status","Sept","Sept Status","Oct","Oct Status","Nov","Nov Status","Dec","Dec Status","Annual","Annual Status","Spring","Spring Status","Summer","Summer Status","Fall","Fall Status","Winter","Winter Status"))

dat2=rbind(dat2,cbind(dat1$Year,dat1$Jan,dat1$Feb,dat1$Mar,dat1$Apr,dat1$May,dat1$Jun,dat1$Jul,dat1$Aug,dat1$Sept,dat1$Oct,dat1$Nov,dat1$Dec,rep(49.78,length(dat1$Year)),rep(-92.83,length(dat1$Year)),rep(413,length(dat1$Year))))

#KENORA
dat1=read.csv("mm6034076.txt",header=FALSE,skip=4,,col.names=c("Year","Jan","Jan Status","Feb","Feb Status","Mar","Mar Status","Apr","Apr Status","May","May Status","Jun","Jun Status","Jul","Jul Status","Aug","Aug Status","Sept","Sept Status","Oct","Oct Status","Nov","Nov Status","Dec","Dec Status","Annual","Annual Status","Spring","Spring Status","Summer","Summer Status","Fall","Fall Status","Winter","Winter Status"))

dat2=rbind(dat2,cbind(dat1$Year,dat1$Jan,dat1$Feb,dat1$Mar,dat1$Apr,dat1$May,dat1$Jun,dat1$Jul,dat1$Aug,dat1$Sept,dat1$Oct,dat1$Nov,dat1$Dec,rep(49.78,length(dat1$Year)),rep(-94.37,length(dat1$Year)),rep(227,length(dat1$Year))))

#SIOUX LOOKOUT
#dat1=read.csv("mm6037776.txt",header=FALSE,skip=4,,col.names=c("Year","Jan","Jan Status","Feb","Feb Status","Mar","Mar Status","Apr","Apr Status","May","May Status","Jun","Jun Status","Jul","Jul Status","Aug","Aug Status","Sept","Sept Status","Oct","Oct Status","Nov","Nov Status","Dec","Dec Status","Annual","Annual Status","Spring","Spring Status","Summer","Summer Status","Fall","Fall Status","Winter","Winter Status"))

#dat2=rbind(dat2,cbind(dat1$Year,dat1$Jan,dat1$Feb,dat1$Mar,dat1$Apr,dat1$May,dat1$Jun,dat1$Jul,dat1$Aug,dat1$Sept,dat1$Oct,dat1$Nov,dat1$Dec,rep(50.12,length(dat1$Year)),rep(-91.9,length(dat1$Year)),rep(383,length(dat1$Year))))

#CAMERON FALLS
dat1=read.csv("mm6041110.txt",header=FALSE,skip=4,,col.names=c("Year","Jan","Jan Status","Feb","Feb Status","Mar","Mar Status","Apr","Apr Status","May","May Status","Jun","Jun Status","Jul","Jul Status","Aug","Aug Status","Sept","Sept Status","Oct","Oct Status","Nov","Nov Status","Dec","Dec Status","Annual","Annual Status","Spring","Spring Status","Summer","Summer Status","Fall","Fall Status","Winter","Winter Status"))

dat2=rbind(dat2,cbind(dat1$Year,dat1$Jan,dat1$Feb,dat1$Mar,dat1$Apr,dat1$May,dat1$Jun,dat1$Jul,dat1$Aug,dat1$Sept,dat1$Oct,dat1$Nov,dat1$Dec,rep(49.15,length(dat1$Year)),rep(-88.35,length(dat1$Year)),rep(233,length(dat1$Year))))

#GERALDTON
dat1=read.csv("mm6042716.txt",header=FALSE,skip=4,,col.names=c("Year","Jan","Jan Status","Feb","Feb Status","Mar","Mar Status","Apr","Apr Status","May","May Status","Jun","Jun Status","Jul","Jul Status","Aug","Aug Status","Sept","Sept Status","Oct","Oct Status","Nov","Nov Status","Dec","Dec Status","Annual","Annual Status","Spring","Spring Status","Summer","Summer Status","Fall","Fall Status","Winter","Winter Status"))

dat2=rbind(dat2,cbind(dat1$Year,dat1$Jan,dat1$Feb,dat1$Mar,dat1$Apr,dat1$May,dat1$Jun,dat1$Jul,dat1$Aug,dat1$Sept,dat1$Oct,dat1$Nov,dat1$Dec,rep(49.78,length(dat1$Year)),rep(-86.93,length(dat1$Year)),rep(349,length(dat1$Year))))

#THUNDER BAY
dat1=read.csv("mm6048260.txt",header=FALSE,skip=4,,col.names=c("Year","Jan","Jan Status","Feb","Feb Status","Mar","Mar Status","Apr","Apr Status","May","May Status","Jun","Jun Status","Jul","Jul Status","Aug","Aug Status","Sept","Sept Status","Oct","Oct Status","Nov","Nov Status","Dec","Dec Status","Annual","Annual Status","Spring","Spring Status","Summer","Summer Status","Fall","Fall Status","Winter","Winter Status"))

dat2=rbind(dat2,cbind(dat1$Year,dat1$Jan,dat1$Feb,dat1$Mar,dat1$Apr,dat1$May,dat1$Jun,dat1$Jul,dat1$Aug,dat1$Sept,dat1$Oct,dat1$Nov,dat1$Dec,rep(48.37,length(dat1$Year)),rep(-89.33,length(dat1$Year)),rep(199,length(dat1$Year))))

#HORNEPAYNE
dat1=read.csv("mm6053575.txt",header=FALSE,skip=4,,col.names=c("Year","Jan","Jan Status","Feb","Feb Status","Mar","Mar Status","Apr","Apr Status","May","May Status","Jun","Jun Status","Jul","Jul Status","Aug","Aug Status","Sept","Sept Status","Oct","Oct Status","Nov","Nov Status","Dec","Dec Status","Annual","Annual Status","Spring","Spring Status","Summer","Summer Status","Fall","Fall Status","Winter","Winter Status"))

dat2=rbind(dat2,cbind(dat1$Year,dat1$Jan,dat1$Feb,dat1$Mar,dat1$Apr,dat1$May,dat1$Jun,dat1$Jul,dat1$Aug,dat1$Sept,dat1$Oct,dat1$Nov,dat1$Dec,rep(49.2,length(dat1$Year)),rep(-84.77,length(dat1$Year)),rep(335,length(dat1$Year))))

#SAULT STE MARIE
dat1=read.csv("mm6057591.txt",header=FALSE,skip=4,,col.names=c("Year","Jan","Jan Status","Feb","Feb Status","Mar","Mar Status","Apr","Apr Status","May","May Status","Jun","Jun Status","Jul","Jul Status","Aug","Aug Status","Sept","Sept Status","Oct","Oct Status","Nov","Nov Status","Dec","Dec Status","Annual","Annual Status","Spring","Spring Status","Summer","Summer Status","Fall","Fall Status","Winter","Winter Status"))

dat2=rbind(dat2,cbind(dat1$Year,dat1$Jan,dat1$Feb,dat1$Mar,dat1$Apr,dat1$May,dat1$Jun,dat1$Jul,dat1$Aug,dat1$Sept,dat1$Oct,dat1$Nov,dat1$Dec,rep(46.48,length(dat1$Year)),rep(-84.52,length(dat1$Year)),rep(192,length(dat1$Year))))

#WAWA
dat1=read.csv("mm6059407.txt",header=FALSE,skip=4,,col.names=c("Year","Jan","Jan Status","Feb","Feb Status","Mar","Mar Status","Apr","Apr Status","May","May Status","Jun","Jun Status","Jul","Jul Status","Aug","Aug Status","Sept","Sept Status","Oct","Oct Status","Nov","Nov Status","Dec","Dec Status","Annual","Annual Status","Spring","Spring Status","Summer","Summer Status","Fall","Fall Status","Winter","Winter Status"))

dat2=rbind(dat2,cbind(dat1$Year,dat1$Jan,dat1$Feb,dat1$Mar,dat1$Apr,dat1$May,dat1$Jun,dat1$Jul,dat1$Aug,dat1$Sept,dat1$Oct,dat1$Nov,dat1$Dec,rep(47.97,length(dat1$Year)),rep(-84.78,length(dat1$Year)),rep(287,length(dat1$Year))))

#CHAPLEAU
dat1=read.csv("mm6061361.txt",header=FALSE,skip=4,,col.names=c("Year","Jan","Jan Status","Feb","Feb Status","Mar","Mar Status","Apr","Apr Status","May","May Status","Jun","Jun Status","Jul","Jul Status","Aug","Aug Status","Sept","Sept Status","Oct","Oct Status","Nov","Nov Status","Dec","Dec Status","Annual","Annual Status","Spring","Spring Status","Summer","Summer Status","Fall","Fall Status","Winter","Winter Status"))

dat2=rbind(dat2,cbind(dat1$Year,dat1$Jan,dat1$Feb,dat1$Mar,dat1$Apr,dat1$May,dat1$Jun,dat1$Jul,dat1$Aug,dat1$Sept,dat1$Oct,dat1$Nov,dat1$Dec,rep(47.82,length(dat1$Year)),rep(-83.35,length(dat1$Year)),rep(447,length(dat1$Year))))

#SUDBURY
dat1=read.csv("mm6068153.txt",header=FALSE,skip=4,,col.names=c("Year","Jan","Jan Status","Feb","Feb Status","Mar","Mar Status","Apr","Apr Status","May","May Status","Jun","Jun Status","Jul","Jul Status","Aug","Aug Status","Sept","Sept Status","Oct","Oct Status","Nov","Nov Status","Dec","Dec Status","Annual","Annual Status","Spring","Spring Status","Summer","Summer Status","Fall","Fall Status","Winter","Winter Status"))

dat2=rbind(dat2,cbind(dat1$Year,dat1$Jan,dat1$Feb,dat1$Mar,dat1$Apr,dat1$May,dat1$Jun,dat1$Jul,dat1$Aug,dat1$Sept,dat1$Oct,dat1$Nov,dat1$Dec,rep(46.62,length(dat1$Year)),rep(-80.8,length(dat1$Year)),rep(348,length(dat1$Year))))

#EARLTON
dat1=read.csv("mm6072230.txt",header=FALSE,skip=4,,col.names=c("Year","Jan","Jan Status","Feb","Feb Status","Mar","Mar Status","Apr","Apr Status","May","May Status","Jun","Jun Status","Jul","Jul Status","Aug","Aug Status","Sept","Sept Status","Oct","Oct Status","Nov","Nov Status","Dec","Dec Status","Annual","Annual Status","Spring","Spring Status","Summer","Summer Status","Fall","Fall Status","Winter","Winter Status"))

dat2=rbind(dat2,cbind(dat1$Year,dat1$Jan,dat1$Feb,dat1$Mar,dat1$Apr,dat1$May,dat1$Jun,dat1$Jul,dat1$Aug,dat1$Sept,dat1$Oct,dat1$Nov,dat1$Dec,rep(47.7,length(dat1$Year)),rep(-79.85,length(dat1$Year)),rep(243,length(dat1$Year))))

#IROQUOIS FALLS
dat1=read.csv("mm6073810.txt",header=FALSE,skip=4,,col.names=c("Year","Jan","Jan Status","Feb","Feb Status","Mar","Mar Status","Apr","Apr Status","May","May Status","Jun","Jun Status","Jul","Jul Status","Aug","Aug Status","Sept","Sept Status","Oct","Oct Status","Nov","Nov Status","Dec","Dec Status","Annual","Annual Status","Spring","Spring Status","Summer","Summer Status","Fall","Fall Status","Winter","Winter Status"))

dat2=rbind(dat2,cbind(dat1$Year,dat1$Jan,dat1$Feb,dat1$Mar,dat1$Apr,dat1$May,dat1$Jun,dat1$Jul,dat1$Aug,dat1$Sept,dat1$Oct,dat1$Nov,dat1$Dec,rep(48.75,length(dat1$Year)),rep(-80.67,length(dat1$Year)),rep(259,length(dat1$Year))))

#KAPUSKASING
dat1=read.csv("mm6073976.txt",header=FALSE,skip=4,,col.names=c("Year","Jan","Jan Status","Feb","Feb Status","Mar","Mar Status","Apr","Apr Status","May","May Status","Jun","Jun Status","Jul","Jul Status","Aug","Aug Status","Sept","Sept Status","Oct","Oct Status","Nov","Nov Status","Dec","Dec Status","Annual","Annual Status","Spring","Spring Status","Summer","Summer Status","Fall","Fall Status","Winter","Winter Status"))

dat2=rbind(dat2,cbind(dat1$Year,dat1$Jan,dat1$Feb,dat1$Mar,dat1$Apr,dat1$May,dat1$Jun,dat1$Jul,dat1$Aug,dat1$Sept,dat1$Oct,dat1$Nov,dat1$Dec,rep(49.42,length(dat1$Year)),rep(-82.47,length(dat1$Year)),rep(227,length(dat1$Year))))

#MOOSONEE
dat1=read.csv("mm6075425.txt",header=FALSE,skip=4,,col.names=c("Year","Jan","Jan Status","Feb","Feb Status","Mar","Mar Status","Apr","Apr Status","May","May Status","Jun","Jun Status","Jul","Jul Status","Aug","Aug Status","Sept","Sept Status","Oct","Oct Status","Nov","Nov Status","Dec","Dec Status","Annual","Annual Status","Spring","Spring Status","Summer","Summer Status","Fall","Fall Status","Winter","Winter Status"))

dat2=rbind(dat2,cbind(dat1$Year,dat1$Jan,dat1$Feb,dat1$Mar,dat1$Apr,dat1$May,dat1$Jun,dat1$Jul,dat1$Aug,dat1$Sept,dat1$Oct,dat1$Nov,dat1$Dec,rep(51.27,length(dat1$Year)),rep(-80.65,length(dat1$Year)),rep(10,length(dat1$Year))))

#SMOKY FALLS
dat1=read.csv("mm6077845.txt",header=FALSE,skip=4,,col.names=c("Year","Jan","Jan Status","Feb","Feb Status","Mar","Mar Status","Apr","Apr Status","May","May Status","Jun","Jun Status","Jul","Jul Status","Aug","Aug Status","Sept","Sept Status","Oct","Oct Status","Nov","Nov Status","Dec","Dec Status","Annual","Annual Status","Spring","Spring Status","Summer","Summer Status","Fall","Fall Status","Winter","Winter Status"))

dat2=rbind(dat2,cbind(dat1$Year,dat1$Jan,dat1$Feb,dat1$Mar,dat1$Apr,dat1$May,dat1$Jun,dat1$Jul,dat1$Aug,dat1$Sept,dat1$Oct,dat1$Nov,dat1$Dec,rep(50.07,length(dat1$Year)),rep(-82.17,length(dat1$Year)),rep(183,length(dat1$Year))))

#TIMMINS
dat1=read.csv("mm6078286.txt",header=FALSE,skip=4,,col.names=c("Year","Jan","Jan Status","Feb","Feb Status","Mar","Mar Status","Apr","Apr Status","May","May Status","Jun","Jun Status","Jul","Jul Status","Aug","Aug Status","Sept","Sept Status","Oct","Oct Status","Nov","Nov Status","Dec","Dec Status","Annual","Annual Status","Spring","Spring Status","Summer","Summer Status","Fall","Fall Status","Winter","Winter Status"))

dat2=rbind(dat2,cbind(dat1$Year,dat1$Jan,dat1$Feb,dat1$Mar,dat1$Apr,dat1$May,dat1$Jun,dat1$Jul,dat1$Aug,dat1$Sept,dat1$Oct,dat1$Nov,dat1$Dec,rep(48.57,length(dat1$Year)),rep(-81.38,length(dat1$Year)),rep(295,length(dat1$Year))))

#MADAWASKA
dat1=read.csv("mm6080192.txt",header=FALSE,skip=4,,col.names=c("Year","Jan","Jan Status","Feb","Feb Status","Mar","Mar Status","Apr","Apr Status","May","May Status","Jun","Jun Status","Jul","Jul Status","Aug","Aug Status","Sept","Sept Status","Oct","Oct Status","Nov","Nov Status","Dec","Dec Status","Annual","Annual Status","Spring","Spring Status","Summer","Summer Status","Fall","Fall Status","Winter","Winter Status"))

dat2=rbind(dat2,cbind(dat1$Year,dat1$Jan,dat1$Feb,dat1$Mar,dat1$Apr,dat1$May,dat1$Jun,dat1$Jul,dat1$Aug,dat1$Sept,dat1$Oct,dat1$Nov,dat1$Dec,rep(45.5,length(dat1$Year)),rep(-77.98,length(dat1$Year)),rep(316,length(dat1$Year))))

#NORTH BAY
dat1=read.csv("mm6085701.txt",header=FALSE,skip=4,,col.names=c("Year","Jan","Jan Status","Feb","Feb Status","Mar","Mar Status","Apr","Apr Status","May","May Status","Jun","Jun Status","Jul","Jul Status","Aug","Aug Status","Sept","Sept Status","Oct","Oct Status","Nov","Nov Status","Dec","Dec Status","Annual","Annual Status","Spring","Spring Status","Summer","Summer Status","Fall","Fall Status","Winter","Winter Status"))

dat2=rbind(dat2,cbind(dat1$Year,dat1$Jan,dat1$Feb,dat1$Mar,dat1$Apr,dat1$May,dat1$Jun,dat1$Jul,dat1$Aug,dat1$Sept,dat1$Oct,dat1$Nov,dat1$Dec,rep(46.37,length(dat1$Year)),rep(-79.42,length(dat1$Year)),rep(370,length(dat1$Year))))

#BROCKVILLE
dat1=read.csv("mm6100971.txt",header=FALSE,skip=4,,col.names=c("Year","Jan","Jan Status","Feb","Feb Status","Mar","Mar Status","Apr","Apr Status","May","May Status","Jun","Jun Status","Jul","Jul Status","Aug","Aug Status","Sept","Sept Status","Oct","Oct Status","Nov","Nov Status","Dec","Dec Status","Annual","Annual Status","Spring","Spring Status","Summer","Summer Status","Fall","Fall Status","Winter","Winter Status"))

dat2=rbind(dat2,cbind(dat1$Year,dat1$Jan,dat1$Feb,dat1$Mar,dat1$Apr,dat1$May,dat1$Jun,dat1$Jul,dat1$Aug,dat1$Sept,dat1$Oct,dat1$Nov,dat1$Dec,rep(44.6,length(dat1$Year)),rep(-75.67,length(dat1$Year)),rep(96,length(dat1$Year))))

#CORNWALL
dat1=read.csv("mm6101874.txt",header=FALSE,skip=4,,col.names=c("Year","Jan","Jan Status","Feb","Feb Status","Mar","Mar Status","Apr","Apr Status","May","May Status","Jun","Jun Status","Jul","Jul Status","Aug","Aug Status","Sept","Sept Status","Oct","Oct Status","Nov","Nov Status","Dec","Dec Status","Annual","Annual Status","Spring","Spring Status","Summer","Summer Status","Fall","Fall Status","Winter","Winter Status"))

dat2=rbind(dat2,cbind(dat1$Year,dat1$Jan,dat1$Feb,dat1$Mar,dat1$Apr,dat1$May,dat1$Jun,dat1$Jul,dat1$Aug,dat1$Sept,dat1$Oct,dat1$Nov,dat1$Dec,rep(45.02,length(dat1$Year)),rep(-74.75,length(dat1$Year)),rep(64,length(dat1$Year))))

#KINGSTON
dat1=read.csv("mm6104142.txt",header=FALSE,skip=4,,col.names=c("Year","Jan","Jan Status","Feb","Feb Status","Mar","Mar Status","Apr","Apr Status","May","May Status","Jun","Jun Status","Jul","Jul Status","Aug","Aug Status","Sept","Sept Status","Oct","Oct Status","Nov","Nov Status","Dec","Dec Status","Annual","Annual Status","Spring","Spring Status","Summer","Summer Status","Fall","Fall Status","Winter","Winter Status"))

dat2=rbind(dat2,cbind(dat1$Year,dat1$Jan,dat1$Feb,dat1$Mar,dat1$Apr,dat1$May,dat1$Jun,dat1$Jul,dat1$Aug,dat1$Sept,dat1$Oct,dat1$Nov,dat1$Dec,rep(44.22,length(dat1$Year)),rep(-76.6,length(dat1$Year)),rep(93,length(dat1$Year))))

#MORRISBURG
dat1=read.csv("mm6105460.txt",header=FALSE,skip=4,,col.names=c("Year","Jan","Jan Status","Feb","Feb Status","Mar","Mar Status","Apr","Apr Status","May","May Status","Jun","Jun Status","Jul","Jul Status","Aug","Aug Status","Sept","Sept Status","Oct","Oct Status","Nov","Nov Status","Dec","Dec Status","Annual","Annual Status","Spring","Spring Status","Summer","Summer Status","Fall","Fall Status","Winter","Winter Status"))

dat2=rbind(dat2,cbind(dat1$Year,dat1$Jan,dat1$Feb,dat1$Mar,dat1$Apr,dat1$May,dat1$Jun,dat1$Jul,dat1$Aug,dat1$Sept,dat1$Oct,dat1$Nov,dat1$Dec,rep(44.92,length(dat1$Year)),rep(-75.18,length(dat1$Year)),rep(82,length(dat1$Year))))

#OTTAWA
dat1=read.csv("mm6105976.txt",header=FALSE,skip=4,,col.names=c("Year","Jan","Jan Status","Feb","Feb Status","Mar","Mar Status","Apr","Apr Status","May","May Status","Jun","Jun Status","Jul","Jul Status","Aug","Aug Status","Sept","Sept Status","Oct","Oct Status","Nov","Nov Status","Dec","Dec Status","Annual","Annual Status","Spring","Spring Status","Summer","Summer Status","Fall","Fall Status","Winter","Winter Status"))

dat2=rbind(dat2,cbind(dat1$Year,dat1$Jan,dat1$Feb,dat1$Mar,dat1$Apr,dat1$May,dat1$Jun,dat1$Jul,dat1$Aug,dat1$Sept,dat1$Oct,dat1$Nov,dat1$Dec,rep(45.38,length(dat1$Year)),rep(-75.72,length(dat1$Year)),rep(79,length(dat1$Year))))

#BEATRICE
#dat1=read.csv("mm6110607.txt",header=FALSE,skip=4,,col.names=c("Year","Jan","Jan Status","Feb","Feb Status","Mar","Mar Status","Apr","Apr Status","May","May Status","Jun","Jun Status","Jul","Jul Status","Aug","Aug Status","Sept","Sept Status","Oct","Oct Status","Nov","Nov Status","Dec","Dec Status","Annual","Annual Status","Spring","Spring Status","Summer","Summer Status","Fall","Fall Status","Winter","Winter Status"))

#dat2=rbind(dat2,cbind(dat1$Year,dat1$Jan,dat1$Feb,dat1$Mar,dat1$Apr,dat1$May,dat1$Jun,dat1$Jul,dat1$Aug,dat1$Sept,dat1$Oct,dat1$Nov,dat1$Dec,rep(45.13,length(dat1$Year)),rep(-79.4,length(dat1$Year)),rep(297,length(dat1$Year))))

#OWEN SOUND
dat1=read.csv("mm6116132.txt",header=FALSE,skip=4,,col.names=c("Year","Jan","Jan Status","Feb","Feb Status","Mar","Mar Status","Apr","Apr Status","May","May Status","Jun","Jun Status","Jul","Jul Status","Aug","Aug Status","Sept","Sept Status","Oct","Oct Status","Nov","Nov Status","Dec","Dec Status","Annual","Annual Status","Spring","Spring Status","Summer","Summer Status","Fall","Fall Status","Winter","Winter Status"))

dat2=rbind(dat2,cbind(dat1$Year,dat1$Jan,dat1$Feb,dat1$Mar,dat1$Apr,dat1$May,dat1$Jun,dat1$Jul,dat1$Aug,dat1$Sept,dat1$Oct,dat1$Nov,dat1$Dec,rep(44.58,length(dat1$Year)),rep(-80.93,length(dat1$Year)),rep(179,length(dat1$Year))))

#HARROW
#dat1=read.csv("mm6130257.txt",header=FALSE,skip=4,,col.names=c("Year","Jan","Jan Status","Feb","Feb Status","Mar","Mar Status","Apr","Apr Status","May","May Status","Jun","Jun Status","Jul","Jul Status","Aug","Aug Status","Sept","Sept Status","Oct","Oct Status","Nov","Nov Status","Dec","Dec Status","Annual","Annual Status","Spring","Spring Status","Summer","Summer Status","Fall","Fall Status","Winter","Winter Status"))

#dat2=rbind(dat2,cbind(dat1$Year,dat1$Jan,dat1$Feb,dat1$Mar,dat1$Apr,dat1$May,dat1$Jun,dat1$Jul,dat1$Aug,dat1$Sept,dat1$Oct,dat1$Nov,dat1$Dec,rep(42.03,length(dat1$Year)),rep(-82.9,length(dat1$Year)),rep(182,length(dat1$Year))))

#RIDGETOWN
dat1=read.csv("mm6137154.txt",header=FALSE,skip=4,,col.names=c("Year","Jan","Jan Status","Feb","Feb Status","Mar","Mar Status","Apr","Apr Status","May","May Status","Jun","Jun Status","Jul","Jul Status","Aug","Aug Status","Sept","Sept Status","Oct","Oct Status","Nov","Nov Status","Dec","Dec Status","Annual","Annual Status","Spring","Spring Status","Summer","Summer Status","Fall","Fall Status","Winter","Winter Status"))

dat2=rbind(dat2,cbind(dat1$Year,dat1$Jan,dat1$Feb,dat1$Mar,dat1$Apr,dat1$May,dat1$Jun,dat1$Jul,dat1$Aug,dat1$Sept,dat1$Oct,dat1$Nov,dat1$Dec,rep(42.45,length(dat1$Year)),rep(-81.88,length(dat1$Year)),rep(206,length(dat1$Year))))

#VINELAND
dat1=read.csv("mm6139148.txt",header=FALSE,skip=4,,col.names=c("Year","Jan","Jan Status","Feb","Feb Status","Mar","Mar Status","Apr","Apr Status","May","May Status","Jun","Jun Status","Jul","Jul Status","Aug","Aug Status","Sept","Sept Status","Oct","Oct Status","Nov","Nov Status","Dec","Dec Status","Annual","Annual Status","Spring","Spring Status","Summer","Summer Status","Fall","Fall Status","Winter","Winter Status"))

dat2=rbind(dat2,cbind(dat1$Year,dat1$Jan,dat1$Feb,dat1$Mar,dat1$Apr,dat1$May,dat1$Jun,dat1$Jul,dat1$Aug,dat1$Sept,dat1$Oct,dat1$Nov,dat1$Dec,rep(43.17,length(dat1$Year)),rep(-79.42,length(dat1$Year)),rep(79,length(dat1$Year))))

#WELLAND
dat1=read.csv("mm6139449.txt",header=FALSE,skip=4,,col.names=c("Year","Jan","Jan Status","Feb","Feb Status","Mar","Mar Status","Apr","Apr Status","May","May Status","Jun","Jun Status","Jul","Jul Status","Aug","Aug Status","Sept","Sept Status","Oct","Oct Status","Nov","Nov Status","Dec","Dec Status","Annual","Annual Status","Spring","Spring Status","Summer","Summer Status","Fall","Fall Status","Winter","Winter Status"))

dat2=rbind(dat2,cbind(dat1$Year,dat1$Jan,dat1$Feb,dat1$Mar,dat1$Apr,dat1$May,dat1$Jun,dat1$Jul,dat1$Aug,dat1$Sept,dat1$Oct,dat1$Nov,dat1$Dec,rep(43,length(dat1$Year)),rep(-79.27,length(dat1$Year)),rep(175,length(dat1$Year))))

#WINDSOR
dat1=read.csv("mm6139527.txt",header=FALSE,skip=4,,col.names=c("Year","Jan","Jan Status","Feb","Feb Status","Mar","Mar Status","Apr","Apr Status","May","May Status","Jun","Jun Status","Jul","Jul Status","Aug","Aug Status","Sept","Sept Status","Oct","Oct Status","Nov","Nov Status","Dec","Dec Status","Annual","Annual Status","Spring","Spring Status","Summer","Summer Status","Fall","Fall Status","Winter","Winter Status"))

dat2=rbind(dat2,cbind(dat1$Year,dat1$Jan,dat1$Feb,dat1$Mar,dat1$Apr,dat1$May,dat1$Jun,dat1$Jul,dat1$Aug,dat1$Sept,dat1$Oct,dat1$Nov,dat1$Dec,rep(42.27,length(dat1$Year)),rep(-82.97,length(dat1$Year)),rep(190,length(dat1$Year))))

#LONDON
dat1=read.csv("mm6144478.txt",header=FALSE,skip=4,,col.names=c("Year","Jan","Jan Status","Feb","Feb Status","Mar","Mar Status","Apr","Apr Status","May","May Status","Jun","Jun Status","Jul","Jul Status","Aug","Aug Status","Sept","Sept Status","Oct","Oct Status","Nov","Nov Status","Dec","Dec Status","Annual","Annual Status","Spring","Spring Status","Summer","Summer Status","Fall","Fall Status","Winter","Winter Status"))

dat2=rbind(dat2,cbind(dat1$Year,dat1$Jan,dat1$Feb,dat1$Mar,dat1$Apr,dat1$May,dat1$Jun,dat1$Jul,dat1$Aug,dat1$Sept,dat1$Oct,dat1$Nov,dat1$Dec,rep(43.03,length(dat1$Year)),rep(-81.15,length(dat1$Year)),rep(278,length(dat1$Year))))

#WOODSTOCK
dat1=read.csv("mm6149625.txt",header=FALSE,skip=4,,col.names=c("Year","Jan","Jan Status","Feb","Feb Status","Mar","Mar Status","Apr","Apr Status","May","May Status","Jun","Jun Status","Jul","Jul Status","Aug","Aug Status","Sept","Sept Status","Oct","Oct Status","Nov","Nov Status","Dec","Dec Status","Annual","Annual Status","Spring","Spring Status","Summer","Summer Status","Fall","Fall Status","Winter","Winter Status"))

dat2=rbind(dat2,cbind(dat1$Year,dat1$Jan,dat1$Feb,dat1$Mar,dat1$Apr,dat1$May,dat1$Jun,dat1$Jul,dat1$Aug,dat1$Sept,dat1$Oct,dat1$Nov,dat1$Dec,rep(43.13,length(dat1$Year)),rep(-80.77,length(dat1$Year)),rep(282,length(dat1$Year))))

#BELLEVILLE
dat1=read.csv("mm6150689.txt",header=FALSE,skip=4,,col.names=c("Year","Jan","Jan Status","Feb","Feb Status","Mar","Mar Status","Apr","Apr Status","May","May Status","Jun","Jun Status","Jul","Jul Status","Aug","Aug Status","Sept","Sept Status","Oct","Oct Status","Nov","Nov Status","Dec","Dec Status","Annual","Annual Status","Spring","Spring Status","Summer","Summer Status","Fall","Fall Status","Winter","Winter Status"))

dat2=rbind(dat2,cbind(dat1$Year,dat1$Jan,dat1$Feb,dat1$Mar,dat1$Apr,dat1$May,dat1$Jun,dat1$Jul,dat1$Aug,dat1$Sept,dat1$Oct,dat1$Nov,dat1$Dec,rep(44.15,length(dat1$Year)),rep(-77.4,length(dat1$Year)),rep(76,length(dat1$Year))))

#HAMILTON
dat1=read.csv("mm6153193.txt",header=FALSE,skip=4,,col.names=c("Year","Jan","Jan Status","Feb","Feb Status","Mar","Mar Status","Apr","Apr Status","May","May Status","Jun","Jun Status","Jul","Jul Status","Aug","Aug Status","Sept","Sept Status","Oct","Oct Status","Nov","Nov Status","Dec","Dec Status","Annual","Annual Status","Spring","Spring Status","Summer","Summer Status","Fall","Fall Status","Winter","Winter Status"))

dat2=rbind(dat2,cbind(dat1$Year,dat1$Jan,dat1$Feb,dat1$Mar,dat1$Apr,dat1$May,dat1$Jun,dat1$Jul,dat1$Aug,dat1$Sept,dat1$Oct,dat1$Nov,dat1$Dec,rep(43.17,length(dat1$Year)),rep(-79.93,length(dat1$Year)),rep(238,length(dat1$Year))))

#ORANGEVILLE
dat1=read.csv("mm6155790.txt",header=FALSE,skip=4,,col.names=c("Year","Jan","Jan Status","Feb","Feb Status","Mar","Mar Status","Apr","Apr Status","May","May Status","Jun","Jun Status","Jul","Jul Status","Aug","Aug Status","Sept","Sept Status","Oct","Oct Status","Nov","Nov Status","Dec","Dec Status","Annual","Annual Status","Spring","Spring Status","Summer","Summer Status","Fall","Fall Status","Winter","Winter Status"))

dat2=rbind(dat2,cbind(dat1$Year,dat1$Jan,dat1$Feb,dat1$Mar,dat1$Apr,dat1$May,dat1$Jun,dat1$Jul,dat1$Aug,dat1$Sept,dat1$Oct,dat1$Nov,dat1$Dec,rep(43.92,length(dat1$Year)),rep(-80.08,length(dat1$Year)),rep(412,length(dat1$Year))))

#TORONTO
dat1=read.csv("mm6158355.txt",header=FALSE,skip=4,,col.names=c("Year","Jan","Jan Status","Feb","Feb Status","Mar","Mar Status","Apr","Apr Status","May","May Status","Jun","Jun Status","Jul","Jul Status","Aug","Aug Status","Sept","Sept Status","Oct","Oct Status","Nov","Nov Status","Dec","Dec Status","Annual","Annual Status","Spring","Spring Status","Summer","Summer Status","Fall","Fall Status","Winter","Winter Status"))

dat2=rbind(dat2,cbind(dat1$Year,dat1$Jan,dat1$Feb,dat1$Mar,dat1$Apr,dat1$May,dat1$Jun,dat1$Jul,dat1$Aug,dat1$Sept,dat1$Oct,dat1$Nov,dat1$Dec,rep(43.67,length(dat1$Year)),rep(-79.4,length(dat1$Year)),rep(113,length(dat1$Year))))

#HALIBURTON
dat1=read.csv("mm6163171.txt",header=FALSE,skip=4,,col.names=c("Year","Jan","Jan Status","Feb","Feb Status","Mar","Mar Status","Apr","Apr Status","May","May Status","Jun","Jun Status","Jul","Jul Status","Aug","Aug Status","Sept","Sept Status","Oct","Oct Status","Nov","Nov Status","Dec","Dec Status","Annual","Annual Status","Spring","Spring Status","Summer","Summer Status","Fall","Fall Status","Winter","Winter Status"))

dat2=rbind(dat2,cbind(dat1$Year,dat1$Jan,dat1$Feb,dat1$Mar,dat1$Apr,dat1$May,dat1$Jun,dat1$Jul,dat1$Aug,dat1$Sept,dat1$Oct,dat1$Nov,dat1$Dec,rep(45.03,length(dat1$Year)),rep(-78.53,length(dat1$Year)),rep(330,length(dat1$Year))))

#PETERBOROUGH
dat1=read.csv("mm6166415.txt",header=FALSE,skip=4,,col.names=c("Year","Jan","Jan Status","Feb","Feb Status","Mar","Mar Status","Apr","Apr Status","May","May Status","Jun","Jun Status","Jul","Jul Status","Aug","Aug Status","Sept","Sept Status","Oct","Oct Status","Nov","Nov Status","Dec","Dec Status","Annual","Annual Status","Spring","Spring Status","Summer","Summer Status","Fall","Fall Status","Winter","Winter Status"))

dat2=rbind(dat2,cbind(dat1$Year,dat1$Jan,dat1$Feb,dat1$Mar,dat1$Apr,dat1$May,dat1$Jun,dat1$Jul,dat1$Aug,dat1$Sept,dat1$Oct,dat1$Nov,dat1$Dec,rep(44.23,length(dat1$Year)),rep(-78.37,length(dat1$Year)),rep(191,length(dat1$Year))))


dat5=dat2[dat2[,1]>=2018,]
dat5=dat5[-c(3,4,12,20,26,28),]

nn=data.frame(x1=rep(2008,5),x2=c(45.88,49.15,51.27,45.38,43.67),x3=c(-82.57,-88.35,-80.65,-75.72,-79.40),x4=c(194,233,10,79,113))
dat5[1,8]=predict(model7,nn)[1]
dat5[7,3]=predict(model2,nn)[2]
dat5[13,13]=predict(model12,nn)[3]
dat5[18,3]=predict(model2,nn)[4]
dat5[25,4]=predict(model2,nn)[5]

#mean prediction
newdata=data.frame(x1=rep(2018,length(dat5[,1])),x2=dat5[,14],x3=dat5[,15],x4=dat5[,16])


newdata1=data.frame(latit=rep(dat5[,14],12),longi=rep(dat5[,15],12),eleva=rep(dat5[,16],12),years=rep(2018,length(dat5[,1])*12),mon=c(rep(1,length(dat5[,1])),rep(2,length(dat5[,1])),rep(3,length(dat5[,1])),rep(4,length(dat5[,1])),rep(5,length(dat5[,1])),rep(6,length(dat5[,1])),rep(7,length(dat5[,1])),rep(8,length(dat5[,1])),rep(9,length(dat5[,1])),rep(10,length(dat5[,1])),rep(11,length(dat5[,1])),rep(12,length(dat5[,1]))),tindex=c((newdata$x1-1978)*12+1,(newdata$x1-1978)*12+2,(newdata$x1-1978)*12+3,(newdata$x1-1978)*12+4,(newdata$x1-1978)*12+5,(newdata$x1-1978)*12+6,(newdata$x1-1978)*12+7,(newdata$x1-1978)*12+8,(newdata$x1-1978)*12+9,(newdata$x1-1978)*12+10,(newdata$x1-1978)*12+11,(newdata$x1-1978)*12+12))


mod3prediction=predict(fit1,newdata1)
mod4prediction=predict(fit2,newdata1)

mod3prediction=matrix(mod3prediction,length(dat5[,1]),length(mod3prediction)/length(dat5[,1]))
mod4prediction=matrix(mod4prediction,length(dat5[,1]),length(mod4prediction)/length(dat5[,1]))

mean1=predict(model1,newdata)
mean2=predict(model2,newdata)
mean3=predict(model3,newdata)
mean4=predict(model4,newdata)
mean5=predict(model5,newdata)
mean6=predict(model6,newdata)
mean7=predict(model7,newdata)
mean8=predict(model8,newdata)
mean9=predict(model9,newdata)
mean10=predict(model10,newdata)
mean11=predict(model11,newdata)
mean12=predict(model12,newdata)

meanmat=cbind(mean1,mean2,mean3,mean4,mean5,mean6,mean7,mean8,mean9,mean10,mean11,mean12)
marginalvar=c(op1$par[6],op2$par[6],op3$par[6],op4$par[6],op5$par[7],op6$par[7],op7$par[7],op8$par[7],op9$par[6],op10$par[6],op11$par[6],op12$par[6])

resipre1<-dat5[,2]-mean1
resipre2<-dat5[,3]-mean2
resipre3<-dat5[,4]-mean3
resipre4<-dat5[,5]-mean4
resipre5<-dat5[,6]-mean5
resipre6<-dat5[,7]-mean6
resipre7<-dat5[,8]-mean7
resipre8<-dat5[,9]-mean8
resipre9<-dat5[,10]-mean9
resipre10<-dat5[,11]-mean10
resipre11<-dat5[,12]-mean11
resipre12<-dat5[,13]-mean12


resipreM=cbind(resipre1,resipre2,resipre3,resipre4,resipre5,resipre6,resipre7,resipre8,resipre9,resipre10,resipre11,resipre12)
month1<-pnorm(resipreM[,1],mean=0,sd =op1$par[6])
month2<-pnorm(resipreM[,2],mean=0,sd =op2$par[6])
month3<-pnorm(resipreM[,3],mean=0,sd =op3$par[6])
month4<-pnorm(resipreM[,4],mean=0,sd =op4$par[6])
month5<-pnorm(resipreM[,5],mean=0,sd =op5$par[7])
month6<-pnorm(resipreM[,6],mean=0,sd =op6$par[7])
month7<-pnorm(resipreM[,7],mean=0,sd =op7$par[7])
month8<-pnorm(resipreM[,8],mean=0,sd =op8$par[7])
month9<-pnorm(resipreM[,9],mean=0,sd =op9$par[6])
month10<-pnorm(resipreM[,10],mean=0,sd =op10$par[6])
month11<-pnorm(resipreM[,11],mean=0,sd =op11$par[6])
month12<-pnorm(resipreM[,12],mean=0,sd =op12$par[6])

p=3

pri<-cbind(month1,month2,month3)

mod1=matrix(0,length(dat5[,1]),12)
mod2=matrix(0,length(dat5[,1]),12)
mod3=matrix(0,length(dat5[,1]),12)
mod4=matrix(0,length(dat5[,1]),12)

percentage1=rep(0,12)
percentage2=rep(0,12)
percentage3=rep(0,12)

for(g in (p+1):12){
  
    RVM5<-RVineMatrix(RVM3$Matrix[(13-g):12,(13-g):12], family = RVM3$family[(13-g):12,(13-g):12], par = RVM3$par[(13-g):12,(13-g):12], par2 = RVM3$par2[(13-g):12,(13-g):12],check.pars = TRUE)
    
    
    RVM4<-RVineMatrix(RVM3$Matrix[(14-g):12,(14-g):12], family = RVM3$family[(14-g):12,(14-g):12], par = RVM3$par[(14-g):12,(14-g):12], par2 = RVM3$par2[(14-g):12,(14-g):12],check.pars = TRUE)
    
    
    dat<-c(rep(0,g))
    
    fitted<-rep(0,length(dat5[,g+1]))
    fitted1<-rep(0,length(dat5[,g+1]))
    fitted2<-rep(0,length(dat5[,g+1]))
    for(j in 1:length(dat5[,g+1])){
      a<-seq(-30,30,0.01)
      b<-rep(0,6001)
      c<-rep(0,6001)
      
      c1<-rep(0,length(dat5[,g+1]))
      c1[j]<-RVinePDF(pri[j,1:(g-1)], RVM4, verbose = TRUE)
      
      
      c2<-function(k){
        cop<-pnorm(k,mean=0,sd =marginalvar[g])
        dat<-c(pri[j,1:(g-1)],cop)
        RVinePDF(dat,RVM5,verbose=TRUE)
      }
      
      f1<-function(k){
        dnorm(k,mean=0,sd =marginalvar[g])
      }
      
      f<-function(k){
        k*f1(k)
      }
      
      
      
      
      #fitted
      for(i in 1:6001){
        b[i]<-c2(a[i])
        c[i]<-f1(a[i])
        d<-a*b*c*0.01/c1[j]
      }
      
      fitted[j]<-sum(d)+meanmat[j,g]
      
      #fitted 1
      
      fitted1[j]<-integrate(f,-Inf,Inf)$value+meanmat[j,g]
      
      
      
      
      
      
      
      
      
      
      
    }
    
    mod1[,g]<-abs(fitted-dat5[,g+1])
    mod2[,g]<-abs(fitted1-dat5[,g+1])
    mod3[,g]<-abs(mod3prediction[,g]-dat5[,g+1])
    mod4[,g]<-abs(mod4prediction[,g]-dat5[,g+1])
    
    pri=cbind(pri,pnorm(fitted-meanmat[,g],mean=0,sd=marginalvar[g]))
    
    count=0
    for(l in 1:length(fitted)){
      if(abs(fitted[l]-dat5[l,g+1])<abs(fitted1[l]-dat5[l,g+1])){
        count=count+1
      }else{
        count=count
      }
    }
    percentage1[g]<-count/length(fitted)
    
    count1=0
    for(i in 1:length(fitted)){
      if(abs(fitted[i]-dat5[i,g+1])<abs(mod3prediction[i,g]-dat5[i,g+1])){
        count1=count1+1
      }else{
        count1=count1
      }
    }
    percentage2[g]<-count1/length(fitted)
    
    count2=0
    for(i in 1:length(fitted)){
      if(abs(fitted[i]-dat5[i,g+1])<abs(mod4prediction[i,g]-dat5[i,g+1])){
        count2=count2+1
      }else{
        count2=count2
      }
    }
    percentage3[g]<-count2/length(fitted)

  
  
  
}


mod1sta=rep(0,27)
mod2sta=rep(0,27)
mod3sta=rep(0,27)
mod4sta=rep(0,27)

percentage1=rep(0,27)
percentage2=rep(0,27)
percentage3=rep(0,27)
for(i in 1:27){
  mod1sta[i]=mean(mod1[i,(p+1):12])
  mod2sta[i]=mean(mod2[i,(p+1):12])
  mod3sta[i]=mean(mod3[i,(p+1):12])
  mod4sta[i]=mean(mod4[i,(p+1):12])
  
  count=0
  count1=0
  count2=0
  for(j in (p+1):12){
    if(mod1[i,j]<=mod2[i,j]){
      count=count+1
    }
    if(mod1[i,j]<=mod3[i,j]){
      count1=count1+1
    }
    if(mod1[i,j]<=mod4[i,j]){
      count2=count2+1
    }
  }
  percentage1[i]=count/10
  percentage2[i]=count1/10
  percentage3[i]=count2/10
}

mod1sta=mod1sta[-c(1,7,18,25)]
mod2sta=mod2sta[-c(1,7,18,25)]
mod3sta=mod3sta[-c(1,7,18,25)]
mod4sta=mod4sta[-c(1,7,18,25)]

predsummary=cbind(mod1sta,mod2sta,mod3sta,mod4sta)
colnames(predsummary)=c("VINE4","MRM","LRM","AR")
boxplot(predsummary,xlab="Model",ylab="MAE")


for(i in 1:length(mod1sta)){
  lines(c(mod1sta[i],mod2sta[i],mod3sta[i],mod4sta[i]),type="b")
}

#line plot for stations
plot(mod3sta,xlab="Stations",type="b",ylim=c(0,5),ylab="MAE",lty=3)
lines(mod1sta,type="b",col="blue",lty=1)
lines(mod2sta,type="b",col="red",lty=2)
lines(mod4sta,type="b",col="green",lty=4)

legend(16, 5, legend=c("VINE4", "MRM","LRM","AR"),
       col=c("blue", "red","black","green"), lty=1:4, cex=0.8)

a=a[-c(1,7,18,25),]
b=b[-c(1,7,18,25),]
mod3new=mod3[-c(1,7,18,25),]
mod4new=mod4[-c(1,7,18,25),]

#line plot for certain stations: first five
plot(a[1,],xlab="Months",type="b",ylab="MAE",lty=1,ylim=c(0,10),col="blue",xaxt="n")
axis(1,at=1:9,lab=seq(4,12,1))
lines(b[1,],type="b",lty=2,col="red")
lines(mod3new[1,4:12],type="b",lty=3,col="black")
lines(mod4new[1,4:12],type="b",lty=4,col="green")

legend(6, 10, legend=c("VINE4", "MRM","LRM","AR"),
       col=c("blue", "red","black","green"), lty=1:4, cex=0.8)

plot(a[2,],xlab="Months",type="b",ylab="MAE",lty=1,ylim=c(0,10),col="blue",xaxt="n")
axis(1,at=1:9,lab=seq(4,12,1))
lines(b[2,],type="b",lty=2,col="red")
lines(mod3new[2,4:12],type="b",lty=3,col="black")
lines(mod4new[2,4:12],type="b",lty=4,col="green")

legend(6, 10, legend=c("VINE4", "MRM","LRM","AR"),
       col=c("blue", "red","black","green"), lty=1:4, cex=0.8)

plot(a[3,],xlab="Months",type="b",ylab="MAE",lty=1,ylim=c(0,10),col="blue",xaxt="n")
axis(1,at=1:9,lab=seq(4,12,1))
lines(b[3,],type="b",lty=2,col="red")
lines(mod3new[3,4:12],type="b",lty=3,col="black")
lines(mod4new[3,4:12],type="b",lty=4,col="green")

legend(6, 10, legend=c("VINE4", "MRM","LRM","AR"),
       col=c("blue", "red","black","green"), lty=1:4, cex=0.8)

plot(a[4,],xlab="Months",type="b",ylab="MAE",lty=1,ylim=c(0,10),col="blue",xaxt="n")
axis(1,at=1:9,lab=seq(4,12,1))
lines(b[4,],type="b",lty=2,col="red")
lines(mod3new[4,4:12],type="b",lty=3,col="black")
lines(mod4new[4,4:12],type="b",lty=4,col="green")

legend(6, 10, legend=c("VINE4", "MRM","LRM","AR"),
       col=c("blue", "red","black","green"), lty=1:4, cex=0.8)


plot(a[7,],xlab="Months",type="b",ylab="MAE",lty=1,ylim=c(0,10),col="blue",xaxt="n")
axis(1,at=1:9,lab=seq(4,12,1))
lines(b[7,],type="b",lty=2,col="red")
lines(mod3new[7,4:12],type="b",lty=3,col="black")
lines(mod4new[7,4:12],type="b",lty=4,col="green")

legend(6, 10, legend=c("VINE4", "MRM","LRM","AR"),
       col=c("blue", "red","black","green"), lty=1:4, cex=0.8)


plot(a[18,],xlab="Months",type="b",ylab="MAE",lty=1,ylim=c(0,10),col="blue",xaxt="n")
axis(1,at=1:9,lab=seq(4,12,1))
lines(b[18,],type="b",lty=2,col="red")
lines(mod3new[18,4:12],type="b",lty=3,col="black")
lines(mod4new[18,4:12],type="b",lty=4,col="green")

legend(6, 10, legend=c("VINE4", "MRM","LRM","AR"),
       col=c("blue", "red","black","green"), lty=1:4, cex=0.8)


plot(a[21,],xlab="Months",type="b",ylab="MAE",lty=1,ylim=c(0,10),col="blue",xaxt="n")
axis(1,at=1:9,lab=seq(4,12,1))
lines(b[21,],type="b",lty=2,col="red")
lines(mod3new[21,4:12],type="b",lty=3,col="black")
lines(mod4new[21,4:12],type="b",lty=4,col="green")

legend(6, 10, legend=c("VINE4", "MRM","LRM","AR"),
       col=c("blue", "red","black","green"), lty=1:4, cex=0.8)
