library(VineCopula)
library(Matrix)
library(matrixcalc)
library(matlib)
library(nlme)
library(SpatioTemporal)

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))))

#ATIKOKAN
#dat1=read.csv("mm6020LPQ.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(-91.58,length(dat1$Year)),rep(442,length(dat1$Year))))

#select the data in 1978-2004
dat3=dat2[dat2[,1]>=1978 & dat2[,1]<=2008,]

#marginal model(Sequential Estimation)
#Jan
y=dat3[dat3[,2]!=-9999.9,2]
x1=dat3[dat3[,2]!=-9999.9,1]
x2=dat3[dat3[,2]!=-9999.9,14]
x3=dat3[dat3[,2]!=-9999.9,15]
x4=dat3[dat3[,2]!=-9999.9,16]
model1=lm(y~x1+x2+x3+x4)

#Feb
y=dat3[dat3[,3]!=-9999.9,3]
x1=dat3[dat3[,3]!=-9999.9,1]
x2=dat3[dat3[,3]!=-9999.9,14]
x3=dat3[dat3[,3]!=-9999.9,15]
x4=dat3[dat3[,3]!=-9999.9,16]
model2=lm(y~x1+x2+x3+x4)

#March
y=dat3[dat3[,4]!=-9999.9,4]
x1=dat3[dat3[,4]!=-9999.9,1]
x2=dat3[dat3[,4]!=-9999.9,14]
x3=dat3[dat3[,4]!=-9999.9,15]
x4=dat3[dat3[,4]!=-9999.9,16]
model3=lm(y~x1+x2+poly(x3,2)+x4)

#April
y=dat3[dat3[,5]!=-9999.9,5]
x1=dat3[dat3[,5]!=-9999.9,1]
x2=dat3[dat3[,5]!=-9999.9,14]
x3=dat3[dat3[,5]!=-9999.9,15]
x4=dat3[dat3[,5]!=-9999.9,16]
model4=lm(y~x1+x2+poly(x3,2)+x4)

#May(linear model not good)
quadlon=poly(dat3[,15],2)

y=dat3[dat3[,6]!=-9999.9,6]
x1=dat3[dat3[,6]!=-9999.9,1]
x2=dat3[dat3[,6]!=-9999.9,14]
x3=dat3[dat3[,6]!=-9999.9,15]
x4=dat3[dat3[,6]!=-9999.9,16]
model5=lm(y~x2+poly(x3,2)+x4)

plot(dat3[dat3[,6]!=-9999.9,1],dat3[dat3[,6]!=-9999.9,6])
plot(dat3[dat3[,6]!=-9999.9,14],dat3[dat3[,6]!=-9999.9,6])
plot(dat3[dat3[,6]!=-9999.9,15],dat3[dat3[,6]!=-9999.9,6])
plot(dat3[dat3[,6]!=-9999.9,16],dat3[dat3[,6]!=-9999.9,6])

#June
y=dat3[dat3[,7]!=-9999.9,7]
x1=dat3[dat3[,7]!=-9999.9,1]
x2=dat3[dat3[,7]!=-9999.9,14]
x3=dat3[dat3[,7]!=-9999.9,15]
x4=dat3[dat3[,7]!=-9999.9,16]
model6=lm(y~x1+x2+poly(x3,2)+x4)

#July
y=dat3[dat3[,8]!=-9999.9,8]
x1=dat3[dat3[,8]!=-9999.9,1]
x2=dat3[dat3[,8]!=-9999.9,14]
x3=dat3[dat3[,8]!=-9999.9,15]
x4=dat3[dat3[,8]!=-9999.9,16]
model7=lm(y~x1+x2+poly(x3,2)+x4)

#Aug
y=dat3[dat3[,9]!=-9999.9,9]
x1=dat3[dat3[,9]!=-9999.9,1]
x2=dat3[dat3[,9]!=-9999.9,14]
x3=dat3[dat3[,9]!=-9999.9,15]
x4=dat3[dat3[,9]!=-9999.9,16]
model8=lm(y~x1+x2+poly(x3,2)+x4)

#Sept
y=dat3[dat3[,10]!=-9999.9,10]
x1=dat3[dat3[,10]!=-9999.9,1]
x2=dat3[dat3[,10]!=-9999.9,14]
x3=dat3[dat3[,10]!=-9999.9,15]
x4=dat3[dat3[,10]!=-9999.9,16]
model9=lm(y~x1+x2+poly(x3,2)+x4)

#Oct
y=dat3[dat3[,11]!=-9999.9,11]
x1=dat3[dat3[,11]!=-9999.9,1]
x2=dat3[dat3[,11]!=-9999.9,14]
x3=dat3[dat3[,11]!=-9999.9,15]
x4=dat3[dat3[,11]!=-9999.9,16]
model10=lm(y~x1+x2+x3+x4)

#Nov
y=dat3[dat3[,12]!=-9999.9,12]
x1=dat3[dat3[,12]!=-9999.9,1]
x2=dat3[dat3[,12]!=-9999.9,14]
x3=dat3[dat3[,12]!=-9999.9,15]
x4=dat3[dat3[,12]!=-9999.9,16]
model11=lm(y~x1+x2+x3+x4)

#Dec
y=dat3[dat3[,13]!=-9999.9,13]
x1=dat3[dat3[,13]!=-9999.9,1]
x2=dat3[dat3[,13]!=-9999.9,14]
x3=dat3[dat3[,13]!=-9999.9,15]
x4=dat3[dat3[,13]!=-9999.9,16]
model12=lm(y~x1+x2+x3+x4)

#missing data imputation
j=2
dat3[dat3[,j]==-9999.9,j]=model1$coefficients[1]+model1$coefficients[2]*dat3[dat3[,j]==-9999.9,1]+model1$coefficients[3]*dat3[dat3[,j]==-9999.9,14]+model1$coefficients[4]*dat3[dat3[,j]==-9999.9,15]+model1$coefficients[5]*dat3[dat3[,j]==-9999.9,16]

j=3


dat3[dat3[,j]==-9999.9,j]=model2$coefficients[1]+model2$coefficients[2]*dat3[dat3[,j]==-9999.9,1]+model2$coefficients[3]*dat3[dat3[,j]==-9999.9,14]+model2$coefficients[4]*dat3[dat3[,j]==-9999.9,15]+model2$coefficients[5]*dat3[dat3[,j]==-9999.9,16]

j=4
dat3[dat3[,j]==-9999.9,j]=model3$coefficients[1]+model3$coefficients[2]*dat3[dat3[,j]==-9999.9,1]+model3$coefficients[3]*dat3[dat3[,j]==-9999.9,14]+model3$coefficients[4]*quadlon[dat3[,j]==-9999.9,1]+model3$coefficients[5]*quadlon[dat3[,j]==-9999.9,2]+model3$coefficients[6]*dat3[dat3[,j]==-9999.9,16]

j=5


dat3[dat3[,j]==-9999.9,j]=model4$coefficients[1]+model4$coefficients[2]*dat3[dat3[,j]==-9999.9,1]+model4$coefficients[3]*dat3[dat3[,j]==-9999.9,14]+model4$coefficients[4]*quadlon[dat3[,j]==-9999.9,1]+model4$coefficients[5]*quadlon[dat3[,j]==-9999.9,2]+model4$coefficients[6]*dat3[dat3[,j]==-9999.9,16]

j=6


dat3[dat3[,j]==-9999.9,j]=model5$coefficients[1]+model5$coefficients[2]*dat3[dat3[,j]==-9999.9,14]+model5$coefficients[3]*quadlon[dat3[,j]==-9999.9,1]+model5$coefficients[4]*quadlon[dat3[,j]==-9999.9,2]+model5$coefficients[5]*dat3[dat3[,j]==-9999.9,16]



j=7
newdata=poly(dat3[dat3[,j]==-9999.9,15],2)
dat3[dat3[,j]==-9999.9,j]=model6$coefficients[1]+model6$coefficients[2]*dat3[dat3[,j]==-9999.9,1]+model6$coefficients[3]*dat3[dat3[,j]==-9999.9,14]+model6$coefficients[4]*quadlon[dat3[,j]==-9999.9,1]+model6$coefficients[5]*quadlon[dat3[,j]==-9999.9,2]+model6$coefficients[6]*dat3[dat3[,j]==-9999.9,16]



j=8
newdata=poly(dat3[dat3[,j]==-9999.9,15],2)


dat3[dat3[,j]==-9999.9,j]=model7$coefficients[1]+model7$coefficients[2]*dat3[dat3[,j]==-9999.9,1]+model7$coefficients[3]*dat3[dat3[,j]==-9999.9,14]+model7$coefficients[4]*quadlon[dat3[,j]==-9999.9,1]+model7$coefficients[5]*quadlon[dat3[,j]==-9999.9,2]+model7$coefficients[6]*dat3[dat3[,j]==-9999.9,16]



j=9
newdata=poly(dat3[dat3[,j]==-9999.9,15],2)
dat3[dat3[,j]==-9999.9,j]=model8$coefficients[1]+model8$coefficients[2]*dat3[dat3[,j]==-9999.9,1]+model8$coefficients[3]*dat3[dat3[,j]==-9999.9,14]+model8$coefficients[4]*quadlon[dat3[,j]==-9999.9,1]+model8$coefficients[5]*quadlon[dat3[,j]==-9999.9,2]+model8$coefficients[6]*dat3[dat3[,j]==-9999.9,16]



j=10
dat3[dat3[,j]==-9999.9,j]=model9$coefficients[1]+model9$coefficients[2]*dat3[dat3[,j]==-9999.9,1]+model9$coefficients[3]*dat3[dat3[,j]==-9999.9,14]+model9$coefficients[4]*quadlon[dat3[,j]==-9999.9,1]+model9$coefficients[5]*quadlon[dat3[,j]==-9999.9,2]+model9$coefficients[6]*dat3[dat3[,j]==-9999.9,16]

j=11
dat3[dat3[,j]==-9999.9,j]=model10$coefficients[1]+model10$coefficients[2]*dat3[dat3[,j]==-9999.9,1]+model10$coefficients[3]*dat3[dat3[,j]==-9999.9,14]+model10$coefficients[4]*dat3[dat3[,j]==-9999.9,15]+model10$coefficients[5]*dat3[dat3[,j]==-9999.9,16]

j=12
dat3[dat3[,j]==-9999.9,j]=model11$coefficients[1]+model11$coefficients[2]*dat3[dat3[,j]==-9999.9,1]+model11$coefficients[3]*dat3[dat3[,j]==-9999.9,14]+model11$coefficients[4]*dat3[dat3[,j]==-9999.9,15]+model11$coefficients[5]*dat3[dat3[,j]==-9999.9,16]

j=13
dat3[dat3[,j]==-9999.9,j]=model12$coefficients[1]+model12$coefficients[2]*dat3[dat3[,j]==-9999.9,1]+model12$coefficients[3]*dat3[dat3[,j]==-9999.9,14]+model12$coefficients[4]*dat3[dat3[,j]==-9999.9,15]+model12$coefficients[5]*dat3[dat3[,j]==-9999.9,16]

#estimation analysis
lik1=function(para){
  loglik=0
  dat4=dat3[dat3[,2]!=-9999.9,2]
  lat=dat3[dat3[,2]!=-9999.9,14]
  lon=dat3[dat3[,2]!=-9999.9,15]
  ele=dat3[dat3[,2]!=-9999.9,16]
  yea=dat3[dat3[,2]!=-9999.9,1]
  for(i in 1:length(dat4)){
    loglik=loglik+dnorm(dat4[i]-para[1]-para[2]*yea[i]-para[3]*lat[i]-para[4]*lon[i]-para[5]*ele[i],mean=0,sd=para[6],log=TRUE)
  }
  -loglik
}

lik2=function(para){
  loglik=0
  dat4=dat3[dat3[,3]!=-9999.9,3]
  lat=dat3[dat3[,3]!=-9999.9,14]
  lon=dat3[dat3[,3]!=-9999.9,15]
  ele=dat3[dat3[,3]!=-9999.9,16]
  yea=dat3[dat3[,3]!=-9999.9,1]
  for(i in 1:length(dat4)){
    loglik=loglik+dnorm(dat4[i]-para[1]-para[2]*yea[i]-para[3]*lat[i]-para[4]*lon[i]-para[5]*ele[i],mean=0,sd=para[6],log=TRUE)
  }
  -loglik
}

lik3=function(para){
  loglik=0
  dat4=dat3[dat3[,4]!=-9999.9,4]
  lat=dat3[dat3[,4]!=-9999.9,14]
  lon=quadlon
  ele=dat3[dat3[,4]!=-9999.9,16]
  yea=dat3[dat3[,4]!=-9999.9,1]
  for(i in 1:length(dat4)){
    loglik=loglik+dnorm(dat4[i]-para[1]-para[2]*yea[i]-para[3]*lat[i]-para[4]*lon[i,1]-para[5]*lon[i,2]-para[6]*ele[i],mean=0,sd=para[7],log=TRUE)
  }
  -loglik
}

lik4=function(para){
  loglik=0
  dat4=dat3[dat3[,5]!=-9999.9,5]
  lat=dat3[dat3[,5]!=-9999.9,14]
  lon=quadlon
  ele=dat3[dat3[,5]!=-9999.9,16]
  yea=dat3[dat3[,5]!=-9999.9,1]
  for(i in 1:length(dat4)){
    loglik=loglik+dnorm(dat4[i]-para[1]-para[2]*yea[i]-para[3]*lat[i]-para[4]*lon[i,1]-para[5]*lon[i,2]-para[6]*ele[i],mean=0,sd=para[7],log=TRUE)
  }
  -loglik
}

lik5=function(para){
  loglik=0
  dat4=dat3[dat3[,6]!=-9999.9,6]
  lat=dat3[dat3[,6]!=-9999.9,14]
  lon=quadlon
  ele=dat3[dat3[,6]!=-9999.9,16]
  yea=dat3[dat3[,6]!=-9999.9,1]
  for(i in 1:length(dat4)){
    loglik=loglik+dnorm(dat4[i]-para[1]-para[2]*lat[i]-para[3]*lon[i,1]-para[4]*lon[i,2]-para[5]*ele[i],mean=0,sd=para[6],log=TRUE)
  }
  -loglik
}

lik6=function(para){
  loglik=0
  dat4=dat3[dat3[,7]!=-9999.9,7]
  lat=dat3[dat3[,7]!=-9999.9,14]
  lon=quadlon
  ele=dat3[dat3[,7]!=-9999.9,16]
  yea=dat3[dat3[,7]!=-9999.9,1]
  for(i in 1:length(dat4)){
    loglik=loglik+dnorm(dat4[i]-para[1]-para[2]*yea[i]-para[3]*lat[i]-para[4]*lon[i,1]-para[5]*lon[i,2]-para[6]*ele[i],mean=0,sd=para[7],log=TRUE)
  }
  -loglik
}

lik7=function(para){
  loglik=0
  dat4=dat3[dat3[,8]!=-9999.9,8]
  lat=dat3[dat3[,8]!=-9999.9,14]
  lon=quadlon
  ele=dat3[dat3[,8]!=-9999.9,16]
  yea=dat3[dat3[,8]!=-9999.9,1]
  for(i in 1:length(dat4)){
    loglik=loglik+dnorm(dat4[i]-para[1]-para[2]*yea[i]-para[3]*lat[i]-para[4]*lon[i,1]-para[5]*lon[i,2]-para[6]*ele[i],mean=0,sd=para[7],log=TRUE)
  }
  -loglik
}

lik8=function(para){
  loglik=0
  dat4=dat3[dat3[,9]!=-9999.9,9]
  lat=dat3[dat3[,9]!=-9999.9,14]
  lon=quadlon
  ele=dat3[dat3[,9]!=-9999.9,16]
  yea=dat3[dat3[,9]!=-9999.9,1]
  for(i in 1:length(dat4)){
    loglik=loglik+dnorm(dat4[i]-para[1]-para[2]*yea[i]-para[3]*lat[i]-para[4]*lon[i,1]-para[5]*lon[i,2]-para[6]*ele[i],mean=0,sd=para[7],log=TRUE)
  }
  -loglik
}

lik9=function(para){
  loglik=0
  dat4=dat3[dat3[,10]!=-9999.9,10]
  lat=dat3[dat3[,10]!=-9999.9,14]
  lon=quadlon
  ele=dat3[dat3[,10]!=-9999.9,16]
  yea=dat3[dat3[,10]!=-9999.9,1]
  for(i in 1:length(dat4)){
    loglik=loglik+dnorm(dat4[i]-para[1]-para[2]*yea[i]-para[3]*lat[i]-para[4]*lon[i,1]-para[5]*lon[i,2]-para[6]*ele[i],mean=0,sd=para[7],log=TRUE)
  }
  -loglik
}

lik10=function(para){
  loglik=0
  dat4=dat3[dat3[,11]!=-9999.9,11]
  lat=dat3[dat3[,11]!=-9999.9,14]
  lon=dat3[dat3[,11]!=-9999.9,15]
  ele=dat3[dat3[,11]!=-9999.9,16]
  yea=dat3[dat3[,11]!=-9999.9,1]
  for(i in 1:length(dat4)){
    loglik=loglik+dnorm(dat4[i]-para[1]-para[2]*yea[i]-para[3]*lat[i]-para[4]*lon[i]-para[5]*ele[i],mean=0,sd=para[6],log=TRUE)
  }
  -loglik
}

lik11=function(para){
  loglik=0
  dat4=dat3[dat3[,12]!=-9999.9,12]
  lat=dat3[dat3[,12]!=-9999.9,14]
  lon=dat3[dat3[,12]!=-9999.9,15]
  ele=dat3[dat3[,12]!=-9999.9,16]
  yea=dat3[dat3[,12]!=-9999.9,1]
  for(i in 1:length(dat4)){
    loglik=loglik+dnorm(dat4[i]-para[1]-para[2]*yea[i]-para[3]*lat[i]-para[4]*lon[i]-para[5]*ele[i],mean=0,sd=para[6],log=TRUE)
  }
  -loglik
}

lik12=function(para){
  loglik=0
  dat4=dat3[dat3[,13]!=-9999.9,13]
  lat=dat3[dat3[,13]!=-9999.9,14]
  lon=dat3[dat3[,13]!=-9999.9,15]
  ele=dat3[dat3[,13]!=-9999.9,16]
  yea=dat3[dat3[,13]!=-9999.9,1]
  for(i in 1:length(dat4)){
    loglik=loglik+dnorm(dat4[i]-para[1]-para[2]*yea[i]-para[3]*lat[i]-para[4]*lon[i]-para[5]*ele[i],mean=0,sd=para[6],log=TRUE)
  }
  -loglik
}

start=c(10,-0.1,0.1,0.1,0.1,2)
lo=c(-Inf,-Inf,-Inf,-Inf,-Inf,0)

start1=c(10,-0.1,0.1,-10,10,0.1,2)
lo1=c(-Inf,-Inf,-Inf,-Inf,-Inf,-Inf,0)

library(nlme)
op1=nlminb(start=start,objective=lik1,lower=lo)
op2=nlminb(start=start,objective=lik2,lower=lo)
op3=nlminb(start=start1,objective=lik3,lower=lo1)
op4=nlminb(start=start1,objective=lik4,lower=lo1)
op5=nlminb(start=start,objective=lik5,lower=lo)
op6=nlminb(start=start1,objective=lik6,lower=lo1)
op7=nlminb(start=start1,objective=lik7,lower=lo1)
op8=nlminb(start=start1,objective=lik8,lower=lo1)
op9=nlminb(start=start1,objective=lik9,lower=lo1)
op10=nlminb(start=start,objective=lik10,lower=lo)
op11=nlminb(start=start,objective=lik11,lower=lo)
op12=nlminb(start=start,objective=lik12,lower=lo)

u1=pnorm(dat3[,2]-op1$par[1]-op1$par[2]*dat3[,1]-op1$par[3]*dat3[,14]-op1$par[4]*dat3[,15]-op1$par[5]*dat3[,16],mean=0,sd=op1$par[6])
u2=pnorm(dat3[,3]-op2$par[1]-op2$par[2]*dat3[,1]-op2$par[3]*dat3[,14]-op2$par[4]*dat3[,15]-op2$par[5]*dat3[,16],mean=0,sd=op2$par[6])
u3=pnorm(dat3[,4]-op3$par[1]-op3$par[2]*dat3[,1]-op3$par[3]*dat3[,14]-op3$par[4]*quadlon[,1]-op3$par[5]*quadlon[,2]-op3$par[6]*dat3[,16],mean=0,sd=op3$par[7])
u4=pnorm(dat3[,5]-op4$par[1]-op4$par[2]*dat3[,1]-op4$par[3]*dat3[,14]-op4$par[4]*quadlon[,1]-op4$par[5]*quadlon[,2]-op4$par[6]*dat3[,16],mean=0,sd=op4$par[7])

u5=pnorm(dat3[,6]-op5$par[1]-op5$par[2]*dat3[,14]-op5$par[3]*quadlon[,1]-op5$par[4]*quadlon[,2]-op5$par[5]*dat3[,16],mean=0,sd=op5$par[6])
u6=pnorm(dat3[,7]-op6$par[1]-op6$par[2]*dat3[,1]-op6$par[3]*dat3[,14]-op6$par[4]*quadlon[,1]-op6$par[5]*quadlon[,2]-op6$par[6]*dat3[,16],mean=0,sd=op6$par[7])
u7=pnorm(dat3[,8]-op7$par[1]-op7$par[2]*dat3[,1]-op7$par[3]*dat3[,14]-op7$par[4]*quadlon[,1]-op7$par[5]*quadlon[,2]-op7$par[6]*dat3[,16],mean=0,sd=op7$par[7])
u8=pnorm(dat3[,9]-op8$par[1]-op8$par[2]*dat3[,1]-op8$par[3]*dat3[,14]-op8$par[4]*quadlon[,1]-op8$par[5]*quadlon[,2]-op8$par[6]*dat3[,16],mean=0,sd=op8$par[7])

u9=pnorm(dat3[,10]-op9$par[1]-op9$par[2]*dat3[,1]-op9$par[3]*dat3[,14]-op9$par[4]*quadlon[,1]-op9$par[5]*quadlon[,2]-op9$par[6]*dat3[,16],mean=0,sd=op9$par[7])
u10=pnorm(dat3[,11]-op10$par[1]-op10$par[2]*dat3[,1]-op10$par[3]*dat3[,14]-op10$par[4]*dat3[,15]-op10$par[5]*dat3[,16],mean=0,sd=op10$par[6])
u11=pnorm(dat3[,12]-op11$par[1]-op11$par[2]*dat3[,1]-op11$par[3]*dat3[,14]-op11$par[4]*dat3[,15]-op11$par[5]*dat3[,16],mean=0,sd=op11$par[6])
u12=pnorm(dat3[,13]-op12$par[1]-op12$par[2]*dat3[,1]-op12$par[3]*dat3[,14]-op12$par[4]*dat3[,15]-op12$par[5]*dat3[,16],mean=0,sd=op12$par[6])

simulate2=cbind(u1,u2,u3,u4,u5,u6,u7,u8,u9,u10,u11,u12)

CVINEMATRIX=c(seq(12,1,-1),0,seq(11,1,-1),0,0,seq(10,1,-1),0,0,0,seq(9,1,-1),0,0,0,0,seq(8,1,-1),rep(0,5),seq(7,1,-1),rep(0,6),seq(6,1,-1),rep(0,7),seq(5,1,-1),rep(0,8),seq(4,1,-1),rep(0,9),seq(3,1,-1),rep(0,10),2,1,rep(0,11),1)
CVINEMATRIX=matrix(CVINEMATRIX,12,12)

RVM3=RVineCopSelect(simulate2,Matrix=CVINEMATRIX,selectioncrit = "BIC",indeptest = TRUE)



likdep=function(para){
  Rvinepar3=matrix(0,12,12)
  co=1
  for(j in 1:12){
    for(i in 1:12){
      if(RVM3$par[i,j]==0){
        Rvinepar3[i,j]=0
      }else{
        Rvinepar3[i,j]=para[co]
        co=co+1
      }
    }
  }
  
  Rvinepar4=matrix(0,12,12)
  co1=0
  for(j in 1:12){
    for(i in 1:12){
      if(RVM3$par2[i,j]==0){
        Rvinepar4[i,j]=0
      }else{
        Rvinepar4[i,j]=para[co1+co]
        co1=co1+1
      }
    }
  }
  
  RVM2<-RVineMatrix(CVINEMATRIX,RVM3$family,Rvinepar3,Rvinepar4,check.pars=FALSE)
  -RVineLogLik(simulate2,RVM2,check.pars = TRUE,calculate.V = FALSE)$loglik
}

#boundary
boun=c(1,-0.99,0.99,0,0,
       2,-0.99,0.99,2.001,30,
       3,0.0001,100,0,0,
       13,0.0001,100,0,0,
       23,-100,-0.0001,0,0,
       33,-100,-0.0001,0,0,
       4,1.0001,100,0,0,
       14,1.0001,100,0,0,
       24,-100,-1.0001,0,0,
       34,-100,-1.0001,0,0,
       5,-100,100,0,0,
       6,1.0001,50,0,0,
       16,1.0001,50,0,0,
       26,-50,-1.0001,0,0,
       36,-50,-1.0001,0,0,
       7,0.0001,7,1.0001,7,
       17,0.0001,7,1.0001,7,
       27,-7,-0.0001,-7,-1.0001,
       37,-7,-0.0001,-7,-1.0001,
       8,1.0001,6,1.0001,8,
       18,1.0001,6,1.0001,8,
       28,-6,-1.0001,-8,-1.0001,
       38,-6,-1.0001,-8,-1.0001,
       9,1.0001,6,0.0001,75,
       19,1.0001,6,0.0001,75,
       29,-6,-1.0001,-75,-0.0001,
       39,-6,-1.0001,-75,-0.0001,
       10,1.0001,8,0.0001,1,
       20,1.0001,8,0.0001,1,
       30,-8,-1.0001,-1,-0.0001,
       40,-8,-1.0001,-1,-0.0001,
       104,1.0001,100,0.0001,0.99,
       114,1.0001,100,0.0001,0.99,
       124,-100,-1.0001,0.0001,0.99,
       134,-100,-1.0001,0.0001,0.99,
       204,1.0001,100,0.0001,0.99,
       214,1.0001,100,0.0001,0.99,
       224,-100,-1.0001,0.0001,0.99,
       234,-100,-1.0001,0.0001,0.99
)

boun=matrix(boun,5,length(boun)/5)
boun=t(boun)


start1=NULL
fam=NULL
lo1=NULL
up1=NULL
for(j in 1:12){
  for(i in 1:12){
    if(RVM3$par[i,j]!=0){
      start1=c(start1,RVM3$par[i,j])
      fam=c(fam,RVM3$family[i,j])
      lo1=c(lo1,boun[boun[,1]==RVM3$family[i,j],2])
      up1=c(up1,boun[boun[,1]==RVM3$family[i,j],3])
    }
  }
}

for(j in 1:12){
  for(i in 1:12){
    if(RVM3$par2[i,j]!=0){
      start1=c(start1,RVM3$par2[i,j])
      lo1=c(lo1,boun[boun[,1]==RVM3$family[i,j],4])
      up1=c(up1,boun[boun[,1]==RVM3$family[i,j],5])
    }
  }
}


out5<-nlminb(objective=likdep,start=start1,control=list(trace=10,eval.max=5000,iter.max=5000),lower=lo1,upper=up1)

lik13=function(para){
  -(lik1(para[1:6])+lik2(para[7:12])+lik3(para[13:19])+lik4(para[20:26])+lik5(para[27:32])+lik6(para[33:39])+lik7(para[40:46])+lik8(para[47:53])+lik9(para[54:60])+lik10(para[61:66])+lik11(para[67:72])+lik12(para[73:78]))
}

out<-fdHess(c(op1$par,op2$par,op3$par,op4$par,op5$par,op6$par,op7$par,op8$par,op9$par,op10$par,op11$par,op12$par),lik13)
H1<-out$Hessian

#number of dependence parameters
ndp=length(as.vector(RVM3$par[RVM3$par!=0]))+length(as.vector(RVM3$par2[RVM3$par2!=0]))
quadlon=poly(dat3[,15],2)


likdep1=function(para){
  Rvinepar3=matrix(0,12,12)
  co=1
  for(j in 1:12){
    for(i in 1:12){
      if(RVM3$par[i,j]==0){
        Rvinepar3[i,j]=0
      }else{
        Rvinepar3[i,j]=para[co]
        co=co+1
      }
    }
  }
  
  Rvinepar4=matrix(0,12,12)
  co1=0
  for(j in 1:12){
    for(i in 1:12){
      if(RVM3$par2[i,j]==0){
        Rvinepar4[i,j]=0
      }else{
        Rvinepar4[i,j]=para[co1+co]
        co1=co1+1
      }
    }
  }
  
  u1<-pnorm(dat3[,2]-para[ndp+1]-para[ndp+2]*dat3[,1]-para[ndp+3]*dat3[,14]-para[ndp+4]*dat3[,15]-para[ndp+5]*dat3[,16],mean=0,sd=para[ndp+6])
  u2<-pnorm(dat3[,3]-para[ndp+7]-para[ndp+8]*dat3[,1]-para[ndp+9]*dat3[,14]-para[ndp+10]*dat3[,15]-para[ndp+11]*dat3[,16],mean=0,sd=para[ndp+12])
  u3<-pnorm(dat3[,4]-para[ndp+13]-para[ndp+14]*dat3[,1]-para[ndp+15]*dat3[,14]-para[ndp+16]*quadlon[,1]-para[ndp+17]*quadlon[,2]-para[ndp+18]*dat3[,16],mean=0,sd=para[ndp+19])
  u4<-pnorm(dat3[,5]-para[ndp+20]-para[ndp+21]*dat3[,1]-para[ndp+22]*dat3[,14]-para[ndp+23]*quadlon[,1]-para[ndp+24]*quadlon[,2]-para[ndp+25]*dat3[,16],mean=0,sd=para[ndp+26])
  u5<-pnorm(dat3[,6]-para[ndp+27]-para[ndp+28]*dat3[,14]-para[ndp+29]*quadlon[,1]-para[ndp+30]*quadlon[,2]-para[ndp+31]*dat3[,16],mean=0,sd=para[ndp+32])
  u6<-pnorm(dat3[,7]-para[ndp+33]-para[ndp+34]*dat3[,1]-para[ndp+35]*dat3[,14]-para[ndp+36]*quadlon[,1]-para[ndp+37]*quadlon[,2]-para[ndp+38]*dat3[,16],mean=0,sd=para[ndp+39])
  u7<-pnorm(dat3[,8]-para[ndp+40]-para[ndp+41]*dat3[,1]-para[ndp+42]*dat3[,14]-para[ndp+43]*quadlon[,1]-para[ndp+44]*quadlon[,2]-para[ndp+45]*dat3[,16],mean=0,sd=para[ndp+46])
  u8<-pnorm(dat3[,9]-para[ndp+47]-para[ndp+48]*dat3[,1]-para[ndp+49]*dat3[,14]-para[ndp+50]*quadlon[,1]-para[ndp+51]*quadlon[,2]-para[ndp+52]*dat3[,16],mean=0,sd=para[ndp+53])
  u9<-pnorm(dat3[,10]-para[ndp+54]-para[ndp+55]*dat3[,1]-para[ndp+56]*dat3[,14]-para[ndp+57]*quadlon[,1]-para[ndp+58]*quadlon[,2]-para[ndp+59]*dat3[,16],mean=0,sd=para[ndp+60])
  u10<-pnorm(dat3[,11]-para[ndp+61]-para[ndp+62]*dat3[,1]-para[ndp+63]*dat3[,14]-para[ndp+64]*dat3[,15]-para[ndp+65]*dat3[,16],mean=0,sd=para[ndp+66])
  u11<-pnorm(dat3[,12]-para[ndp+67]-para[ndp+68]*dat3[,1]-para[ndp+69]*dat3[,14]-para[ndp+70]*dat3[,15]-para[ndp+71]*dat3[,16],mean=0,sd=para[ndp+72])
  u12<-pnorm(dat3[,13]-para[ndp+73]-para[ndp+74]*dat3[,1]-para[ndp+75]*dat3[,14]-para[ndp+76]*dat3[,15]-para[ndp+77]*dat3[,16],mean=0,sd=para[ndp+78])
  simulate5<-cbind(u1,u2,u3,u4,u5,u6,u7,u8,u9,u10,u11,u12)
  
  RVM2<-RVineMatrix(CVINEMATRIX,RVM3$family,Rvinepar3,Rvinepar4,check.pars=FALSE)
  RVineLogLik(simulate5,RVM2,check.pars = FALSE,calculate.V = FALSE)$loglik
}

out1<-fdHess(c(out5$par,op1$par,op2$par,op3$par,op4$par,op5$par,op6$par,op7$par,op8$par,op9$par,op10$par,op11$par,op12$par),likdep1)
H2<-out1$Hessian

Q<-H2[1:ndp,(ndp+1):(ndp+78)]%*%inv(H1)


J<-matrix(0,78,78)
omega<-matrix(0,ndp,ndp)

for(i in 1:length(dat3[,1])){
  lik14=function(para){
    dnorm(dat3[i,2]-para[1]-para[2]*dat3[i,1]-para[3]*dat3[i,14]-para[4]*dat3[i,15]-para[5]*dat3[i,16],mean=0,sd=para[6],log=TRUE)+dnorm(dat3[i,3]-para[7]-para[8]*dat3[i,1]-para[9]*dat3[i,14]-para[10]*dat3[i,15]-para[11]*dat3[i,16],mean=0,sd=para[12],log=TRUE)+dnorm(dat3[i,4]-para[13]-para[14]*dat3[i,1]-para[15]*dat3[i,14]--para[16]*quadlon[i,1]-para[17]*quadlon[i,2]-para[18]*dat3[i,16],mean=0,sd=para[19],log=TRUE)+dnorm(dat3[i,5]-para[20]-para[21]*dat3[i,1]-para[22]*dat3[i,14]--para[23]*quadlon[i,1]-para[24]*quadlon[i,2]-para[25]*dat3[i,16],mean=0,sd=para[26],log=TRUE) +dnorm(dat3[i,6]-para[27]-para[28]*dat3[i,14]-para[29]*quadlon[i,1]-para[30]*quadlon[i,2]-para[31]*dat3[i,16],mean=0,sd=para[32],log=TRUE)+dnorm(dat3[i,7]-para[33]-para[34]*dat3[i,1]-para[35]*dat3[i,14]-para[36]*quadlon[i,1]-para[37]*quadlon[i,2]-para[38]*dat3[i,16],mean=0,sd=para[39],log=TRUE)+dnorm(dat3[i,8]-para[40]-para[41]*dat3[i,1]-para[42]*dat3[i,14]-para[43]*quadlon[i,1]-para[44]*quadlon[i,2]-para[45]*dat3[i,16],mean=0,sd=para[46],log=TRUE)+dnorm(dat3[i,9]-para[47]-para[48]*dat3[i,1]-para[49]*dat3[i,14]-para[50]*quadlon[i,1]-para[51]*quadlon[i,2]-para[52]*dat3[i,16],mean=0,sd=para[53],log=TRUE)+dnorm(dat3[i,10]-para[54]-para[55]*dat3[i,1]-para[56]*dat3[i,14]--para[57]*quadlon[i,1]-para[58]*quadlon[i,2]-para[59]*dat3[i,16],mean=0,sd=para[60],log=TRUE)+dnorm(dat3[i,11]-para[61]-para[62]*dat3[i,1]-para[63]*dat3[i,14]-para[64]*dat3[i,15]-para[65]*dat3[i,16],mean=0,sd=para[66],log=TRUE)+dnorm(dat3[i,12]-para[67]-para[68]*dat3[i,1]-para[69]*dat3[i,14]-para[70]*dat3[i,15]-para[71]*dat3[i,16],mean=0,sd=para[72],log=TRUE)+dnorm(dat3[i,13]-para[73]-para[74]*dat3[i,1]-para[75]*dat3[i,14]-para[76]*dat3[i,15]-para[77]*dat3[i,16],mean=0,sd=para[78],log=TRUE)
  }
  out6<-genGradient(c(op1$par,op2$par,op3$par,op4$par,op5$par,op6$par,op7$par,op8$par,op9$par,op10$par,op11$par,op12$par),lik14)
  J<-J+out6%*%t(out6)
  
  lik15=function(para){
    Rvinepar3=matrix(0,12,12)
    co=1
    for(j in 1:12){
      for(k in 1:12){
        if(RVM3$par[k,j]==0){
          Rvinepar3[k,j]=0
        }else{
          Rvinepar3[k,j]=para[co]
          co=co+1
        }
      }
    }
    
    Rvinepar4=matrix(0,12,12)
    co1=0
    for(j in 1:12){
      for(k in 1:12){
        if(RVM3$par2[k,j]==0){
          Rvinepar4[k,j]=0
        }else{
          Rvinepar4[k,j]=para[co1+co]
          co1=co1+1
        }
      }
    }
    
    u1<-pnorm(dat3[i,2]-para[ndp+1]-para[ndp+2]*dat3[i,1]-para[ndp+3]*dat3[i,14]-para[ndp+4]*dat3[i,15]-para[ndp+5]*dat3[i,16],mean=0,sd=para[ndp+6])
    u2<-pnorm(dat3[i,3]-para[ndp+7]-para[ndp+8]*dat3[i,1]-para[ndp+9]*dat3[i,14]-para[ndp+10]*dat3[i,15]-para[ndp+11]*dat3[i,16],mean=0,sd=para[ndp+12])
    u3<-pnorm(dat3[i,4]-para[ndp+13]-para[ndp+14]*dat3[i,1]-para[ndp+15]*dat3[i,14]-para[ndp+16]*quadlon[i,1]-para[ndp+17]*quadlon[i,2]-para[ndp+18]*dat3[i,16],mean=0,sd=para[ndp+19])
    u4<-pnorm(dat3[i,5]-para[ndp+20]-para[ndp+21]*dat3[i,1]-para[ndp+22]*dat3[i,14]-para[ndp+23]*quadlon[i,1]-para[ndp+24]*quadlon[i,2]-para[ndp+25]*dat3[i,16],mean=0,sd=para[ndp+26])
    u5<-pnorm(dat3[i,6]-para[ndp+27]-para[ndp+28]*dat3[i,14]-para[ndp+29]*quadlon[i,1]-para[ndp+30]*quadlon[i,2]-para[ndp+31]*dat3[i,16],mean=0,sd=para[ndp+32])
    u6<-pnorm(dat3[i,7]-para[ndp+33]-para[ndp+34]*dat3[i,1]-para[ndp+35]*dat3[i,14]-para[ndp+36]*quadlon[i,1]-para[ndp+37]*quadlon[i,2]-para[ndp+38]*dat3[i,16],mean=0,sd=para[ndp+39])
    u7<-pnorm(dat3[i,8]-para[ndp+40]-para[ndp+41]*dat3[i,1]-para[ndp+42]*dat3[i,14]-para[ndp+43]*quadlon[i,1]-para[ndp+44]*quadlon[i,2]-para[ndp+45]*dat3[i,16],mean=0,sd=para[ndp+46])
    u8<-pnorm(dat3[i,9]-para[ndp+47]-para[ndp+48]*dat3[i,1]-para[ndp+49]*dat3[i,14]-para[ndp+50]*quadlon[i,1]-para[ndp+51]*quadlon[i,2]-para[ndp+52]*dat3[i,16],mean=0,sd=para[ndp+53])
    u9<-pnorm(dat3[i,10]-para[ndp+54]-para[ndp+55]*dat3[i,1]-para[ndp+56]*dat3[i,14]-para[ndp+57]*quadlon[i,1]-para[ndp+58]*quadlon[i,2]-para[ndp+59]*dat3[i,16],mean=0,sd=para[ndp+60])
    u10<-pnorm(dat3[i,11]-para[ndp+61]-para[ndp+62]*dat3[i,1]-para[ndp+63]*dat3[i,14]-para[ndp+64]*dat3[i,15]-para[ndp+65]*dat3[i,16],mean=0,sd=para[ndp+66])
    u11<-pnorm(dat3[i,12]-para[ndp+67]-para[ndp+68]*dat3[i,1]-para[ndp+69]*dat3[i,14]-para[ndp+70]*dat3[i,15]-para[ndp+71]*dat3[i,16],mean=0,sd=para[ndp+72])
    u12<-pnorm(dat3[i,13]-para[ndp+73]-para[ndp+74]*dat3[i,1]-para[ndp+75]*dat3[i,14]-para[ndp+76]*dat3[i,15]-para[ndp+77]*dat3[i,16],mean=0,sd=para[ndp+78])
    simulate4<-cbind(u1,u2,u3,u4,u5,u6,u7,u8,u9,u10,u11,u12)
    
    RVM2<-RVineMatrix(CVINEMATRIX,RVM3$family,Rvinepar3,Rvinepar4,check.pars=FALSE)
    RVineLogLik(simulate4,RVM2,check.pars = FALSE,calculate.V = FALSE)$loglik
    
  }
  out7<-genGradient(c(out5$par,op1$par,op2$par,op3$par,op4$par,op5$par,op6$par,op7$par,op8$par,op9$par,op10$par,op11$par,op12$par),lik15)
  ome<-out7[1:ndp]-Q%*%out6
  omega<-omega+ome%*%t(ome)
  
}


J<-J/length(dat3[,1])
H1<-H1/length(dat3[,1])
robustase<-inv(H1)%*%J%*%inv(H1)

omega<-omega/length(dat3[,1])
H3<-H2/length(dat3[,1])
robustase1<-inv(H3[1:ndp,1:ndp])%*%omega%*%inv(H3[1:ndp,1:ndp])

(sqrt(diag(robustase))/sqrt(length(dat3[,1])))
(sqrt(diag(robustase1))/sqrt(length(dat3[,1])))










#simultaneous estimation


lik=function(para){
  Rvinepar3=matrix(0,12,12)
  co=1
  for(j in 1:12){
    for(i in 1:12){
      if(RVM3$par[i,j]==0){
        Rvinepar3[i,j]=0
      }else{
        Rvinepar3[i,j]=para[co]
        co=co+1
      }
    }
  }
  
  Rvinepar4=matrix(0,12,12)
  co1=0
  for(j in 1:12){
    for(i in 1:12){
      if(RVM3$par2[i,j]==0){
        Rvinepar4[i,j]=0
      }else{
        Rvinepar4[i,j]=para[co1+co]
        co1=co1+1
      }
    }
  }
  
  u1<-pnorm(dat3[,2]-para[ndp+1]-para[ndp+2]*dat3[,1]-para[ndp+3]*dat3[,14]-para[ndp+4]*dat3[,15]-para[ndp+5]*dat3[,16],mean=0,sd=para[ndp+6])
  u2<-pnorm(dat3[,3]-para[ndp+7]-para[ndp+8]*dat3[,1]-para[ndp+9]*dat3[,14]-para[ndp+10]*dat3[,15]-para[ndp+11]*dat3[,16],mean=0,sd=para[ndp+12])
  u3<-pnorm(dat3[,4]-para[ndp+13]-para[ndp+14]*dat3[,1]-para[ndp+15]*dat3[,14]-para[ndp+16]*dat3[,15]-para[ndp+17]*dat3[,16],mean=0,sd=para[ndp+18])
  u4<-pnorm(dat3[,5]-para[ndp+19]-para[ndp+20]*dat3[,1]-para[ndp+21]*dat3[,14]-para[ndp+22]*dat3[,15]-para[ndp+23]*dat3[,16],mean=0,sd=para[ndp+24])
  u5<-pnorm(dat3[,6]-para[ndp+25]-para[ndp+26]*dat3[,1]-para[ndp+27]*dat3[,14]-para[ndp+28]*quadlon[,1]-para[ndp+29]*quadlon[,2]-para[ndp+30]*dat3[,16],mean=0,sd=para[ndp+31])
  u6<-pnorm(dat3[,7]-para[ndp+32]-para[ndp+33]*dat3[,1]-para[ndp+34]*dat3[,14]-para[ndp+35]*quadlon[,1]-para[ndp+36]*quadlon[,2]-para[ndp+37]*dat3[,16],mean=0,sd=para[ndp+38])
  u7<-pnorm(dat3[,8]-para[ndp+39]-para[ndp+40]*dat3[,1]-para[ndp+41]*dat3[,14]-para[ndp+42]*quadlon[,1]-para[ndp+43]*quadlon[,2]-para[ndp+44]*dat3[,16],mean=0,sd=para[ndp+45])
  u8<-pnorm(dat3[,9]-para[ndp+46]-para[ndp+47]*dat3[,1]-para[ndp+48]*dat3[,14]-para[ndp+49]*quadlon[,1]-para[ndp+50]*quadlon[,2]-para[ndp+51]*dat3[,16],mean=0,sd=para[ndp+52])
  u9<-pnorm(dat3[,10]-para[ndp+53]-para[ndp+54]*dat3[,1]-para[ndp+55]*dat3[,14]-para[ndp+56]*dat3[,15]-para[ndp+57]*dat3[,16],mean=0,sd=para[ndp+58])
  u10<-pnorm(dat3[,11]-para[ndp+59]-para[ndp+60]*dat3[,1]-para[ndp+61]*dat3[,14]-para[ndp+62]*dat3[,15]-para[ndp+63]*dat3[,16],mean=0,sd=para[ndp+64])
  u11<-pnorm(dat3[,12]-para[ndp+65]-para[ndp+66]*dat3[,1]-para[ndp+67]*dat3[,14]-para[ndp+68]*dat3[,15]-para[ndp+69]*dat3[,16],mean=0,sd=para[ndp+70])
  u12<-pnorm(dat3[,13]-para[ndp+71]-para[ndp+72]*dat3[,1]-para[ndp+73]*dat3[,14]-para[ndp+74]*dat3[,15]-para[ndp+75]*dat3[,16],mean=0,sd=para[ndp+76])
  simulate3<-cbind(u1,u2,u3,u4,u5,u6,u7,u8,u9,u10,u11,u12)
  
  RVM2<-RVineMatrix(CVINEMATRIX,RVM3$family,Rvinepar3,Rvinepar4,check.pars=FALSE)
  -(RVineLogLik(simulate3,RVM2,check.pars = TRUE,calculate.V = FALSE)$loglik+sum(dnorm(dat3[,2]-para[ndp+1]-para[ndp+2]*dat3[,1]-para[ndp+3]*dat3[,14]-para[ndp+4]*dat3[,15]-para[ndp+5]*dat3[,16],mean=0,sd=para[ndp+6],log=TRUE))+sum(dnorm(dat3[,3]-para[ndp+7]-para[ndp+8]*dat3[,1]-para[ndp+9]*dat3[,14]-para[ndp+10]*dat3[,15]-para[ndp+11]*dat3[,16],mean=0,sd=para[ndp+12],log=TRUE))+sum(dnorm(dat3[,4]-para[ndp+13]-para[ndp+14]*dat3[,1]-para[ndp+15]*dat3[,14]-para[ndp+16]*dat3[,15]-para[ndp+17]*dat3[,16],mean=0,sd=para[ndp+18],log=TRUE))+sum(dnorm(dat3[,5]-para[ndp+19]-para[ndp+20]*dat3[,1]-para[ndp+21]*dat3[,14]-para[ndp+22]*dat3[,15]-para[ndp+23]*dat3[,16],mean=0,sd=para[ndp+24],log=TRUE)) +sum(dnorm(dat3[,6]-para[ndp+25]-para[ndp+26]*dat3[,1]-para[ndp+27]*dat3[,14]-para[ndp+28]*quadlon[,1]-para[ndp+29]*quadlon[,2]-para[ndp+30]*dat3[,16],mean=0,sd=para[ndp+31],log=TRUE))+sum(dnorm(dat3[,7]-para[ndp+32]-para[ndp+33]*dat3[,1]-para[ndp+34]*dat3[,14]-para[ndp+35]*quadlon[,1]-para[ndp+36]*quadlon[,2]-para[ndp+37]*dat3[,16],mean=0,sd=para[ndp+38],log=TRUE))+sum(dnorm(dat3[,8]-para[ndp+39]-para[ndp+40]*dat3[,1]-para[ndp+41]*dat3[,14]-para[ndp+42]*quadlon[,1]-para[ndp+43]*quadlon[,2]-para[ndp+44]*dat3[,16],mean=0,sd=para[ndp+45],log=TRUE))+sum(dnorm(dat3[,9]-para[ndp+46]-para[ndp+47]*dat3[,1]-para[ndp+48]*dat3[,14]-para[ndp+49]*quadlon[,1]-para[ndp+50]*quadlon[,2]-para[ndp+51]*dat3[,16],mean=0,sd=para[ndp+52],log=TRUE))+sum(dnorm(dat3[,10]-para[ndp+53]-para[ndp+54]*dat3[,1]-para[ndp+55]*dat3[,14]-para[ndp+56]*dat3[,15]-para[ndp+57]*dat3[,16],mean=0,sd=para[ndp+58],log=TRUE))+sum(dnorm(dat3[,11]-para[ndp+59]-para[ndp+60]*dat3[,1]-para[ndp+61]*dat3[,14]-para[ndp+62]*dat3[,15]-para[ndp+63]*dat3[,16],mean=0,sd=para[ndp+64],log=TRUE))+sum(dnorm(dat3[,12]-para[ndp+65]-para[ndp+66]*dat3[,1]-para[ndp+67]*dat3[,14]-para[ndp+68]*dat3[,15]-para[ndp+69]*dat3[,16],mean=0,sd=para[ndp+70],log=TRUE))+sum(dnorm(dat3[,13]-para[ndp+71]-para[ndp+72]*dat3[,1]-para[ndp+73]*dat3[,14]-para[ndp+74]*dat3[,15]-para[ndp+75]*dat3[,16],mean=0,sd=para[ndp+76],log=TRUE)))
  
}

start2=c(start1,op1$par,op2$par,op3$par,op4$par,op5$par,op6$par,op7$par,op8$par,op9$par,op10$par,op11$par,op12$par)
lo3=c(-Inf,-Inf,-Inf,-Inf,-Inf,-Inf,0)
up3=c(Inf,Inf,Inf,Inf,Inf,Inf,Inf)

lo4=c(-Inf,-Inf,-Inf,-Inf,-Inf,0)
up4=c(Inf,Inf,Inf,Inf,Inf,Inf)

lo2=c(lo1,rep(lo4,4),rep(lo3,4),rep(lo4,4))
up2=c(up1,rep(up4,4),rep(up3,4),rep(up4,4))

out111=nlminb(objective=lik,start=start2,control=list(trace=10,eval.max=5000,iter.max=5000),lower=lo2,upper=up2)
out113=nlminb(objective=lik,start=out111$par,control=list(trace=10,eval.max=5000,iter.max=5000),lower=lo2,upper=up2)

likhess=function(para){
  Rvinepar3=matrix(0,12,12)
  co=1
  for(j in 1:12){
    for(i in 1:12){
      if(RVM3$par[i,j]==0){
        Rvinepar3[i,j]=0
      }else{
        Rvinepar3[i,j]=para[co]
        co=co+1
      }
    }
  }
  
  Rvinepar4=matrix(0,12,12)
  co1=0
  for(j in 1:12){
    for(i in 1:12){
      if(RVM3$par2[i,j]==0){
        Rvinepar4[i,j]=0
      }else{
        Rvinepar4[i,j]=para[co1+co]
        co1=co1+1
      }
    }
  }
  
  u1<-pnorm(dat3[,2]-para[ndp+1]-para[ndp+2]*dat3[,1]-para[ndp+3]*dat3[,14]-para[ndp+4]*dat3[,15]-para[ndp+5]*dat3[,16],mean=0,sd=para[ndp+6])
  u2<-pnorm(dat3[,3]-para[ndp+7]-para[ndp+8]*dat3[,1]-para[ndp+9]*dat3[,14]-para[ndp+10]*dat3[,15]-para[ndp+11]*dat3[,16],mean=0,sd=para[ndp+12])
  u3<-pnorm(dat3[,4]-para[ndp+13]-para[ndp+14]*dat3[,1]-para[ndp+15]*dat3[,14]-para[ndp+16]*dat3[,15]-para[ndp+17]*dat3[,16],mean=0,sd=para[ndp+18])
  u4<-pnorm(dat3[,5]-para[ndp+19]-para[ndp+20]*dat3[,1]-para[ndp+21]*dat3[,14]-para[ndp+22]*dat3[,15]-para[ndp+23]*dat3[,16],mean=0,sd=para[ndp+24])
  u5<-pnorm(dat3[,6]-para[ndp+25]-para[ndp+26]*dat3[,1]-para[ndp+27]*dat3[,14]-para[ndp+28]*quadlon[,1]-para[ndp+29]*quadlon[,2]-para[ndp+30]*dat3[,16],mean=0,sd=para[ndp+31])
  u6<-pnorm(dat3[,7]-para[ndp+32]-para[ndp+33]*dat3[,1]-para[ndp+34]*dat3[,14]-para[ndp+35]*quadlon[,1]-para[ndp+36]*quadlon[,2]-para[ndp+37]*dat3[,16],mean=0,sd=para[ndp+38])
  u7<-pnorm(dat3[,8]-para[ndp+39]-para[ndp+40]*dat3[,1]-para[ndp+41]*dat3[,14]-para[ndp+42]*quadlon[,1]-para[ndp+43]*quadlon[,2]-para[ndp+44]*dat3[,16],mean=0,sd=para[ndp+45])
  u8<-pnorm(dat3[,9]-para[ndp+46]-para[ndp+47]*dat3[,1]-para[ndp+48]*dat3[,14]-para[ndp+49]*quadlon[,1]-para[ndp+50]*quadlon[,2]-para[ndp+51]*dat3[,16],mean=0,sd=para[ndp+52])
  u9<-pnorm(dat3[,10]-para[ndp+53]-para[ndp+54]*dat3[,1]-para[ndp+55]*dat3[,14]-para[ndp+56]*dat3[,15]-para[ndp+57]*dat3[,16],mean=0,sd=para[ndp+58])
  u10<-pnorm(dat3[,11]-para[ndp+59]-para[ndp+60]*dat3[,1]-para[ndp+61]*dat3[,14]-para[ndp+62]*dat3[,15]-para[ndp+63]*dat3[,16],mean=0,sd=para[ndp+64])
  u11<-pnorm(dat3[,12]-para[ndp+65]-para[ndp+66]*dat3[,1]-para[ndp+67]*dat3[,14]-para[ndp+68]*dat3[,15]-para[ndp+69]*dat3[,16],mean=0,sd=para[ndp+70])
  u12<-pnorm(dat3[,13]-para[ndp+71]-para[ndp+72]*dat3[,1]-para[ndp+73]*dat3[,14]-para[ndp+74]*dat3[,15]-para[ndp+75]*dat3[,16],mean=0,sd=para[ndp+76])
  simulate3<-cbind(u1,u2,u3,u4,u5,u6,u7,u8,u9,u10,u11,u12)
  
  RVM2<-RVineMatrix(CVINEMATRIX,RVM3$family,Rvinepar3,Rvinepar4,check.pars=FALSE)
  RVineLogLik(simulate3,RVM2,check.pars = TRUE,calculate.V = FALSE)$loglik+sum(dnorm(dat3[,2]-para[ndp+1]-para[ndp+2]*dat3[,1]-para[ndp+3]*dat3[,14]-para[ndp+4]*dat3[,15]-para[ndp+5]*dat3[,16],mean=0,sd=para[ndp+6],log=TRUE))+sum(dnorm(dat3[,3]-para[ndp+7]-para[ndp+8]*dat3[,1]-para[ndp+9]*dat3[,14]-para[ndp+10]*dat3[,15]-para[ndp+11]*dat3[,16],mean=0,sd=para[ndp+12],log=TRUE))+sum(dnorm(dat3[,4]-para[ndp+13]-para[ndp+14]*dat3[,1]-para[ndp+15]*dat3[,14]-para[ndp+16]*dat3[,15]-para[ndp+17]*dat3[,16],mean=0,sd=para[ndp+18],log=TRUE))+sum(dnorm(dat3[,5]-para[ndp+19]-para[ndp+20]*dat3[,1]-para[ndp+21]*dat3[,14]-para[ndp+22]*dat3[,15]-para[ndp+23]*dat3[,16],mean=0,sd=para[ndp+24],log=TRUE)) +sum(dnorm(dat3[,6]-para[ndp+25]-para[ndp+26]*dat3[,1]-para[ndp+27]*dat3[,14]-para[ndp+28]*quadlon[,1]-para[ndp+29]*quadlon[,2]-para[ndp+30]*dat3[,16],mean=0,sd=para[ndp+31],log=TRUE))+sum(dnorm(dat3[,7]-para[ndp+32]-para[ndp+33]*dat3[,1]-para[ndp+34]*dat3[,14]-para[ndp+35]*quadlon[,1]-para[ndp+36]*quadlon[,2]-para[ndp+37]*dat3[,16],mean=0,sd=para[ndp+38],log=TRUE))+sum(dnorm(dat3[,8]-para[ndp+39]-para[ndp+40]*dat3[,1]-para[ndp+41]*dat3[,14]-para[ndp+42]*quadlon[,1]-para[ndp+43]*quadlon[,2]-para[ndp+44]*dat3[,16],mean=0,sd=para[ndp+45],log=TRUE))+sum(dnorm(dat3[,9]-para[ndp+46]-para[ndp+47]*dat3[,1]-para[ndp+48]*dat3[,14]-para[ndp+49]*quadlon[,1]-para[ndp+50]*quadlon[,2]-para[ndp+51]*dat3[,16],mean=0,sd=para[ndp+52],log=TRUE))+sum(dnorm(dat3[,10]-para[ndp+53]-para[ndp+54]*dat3[,1]-para[ndp+55]*dat3[,14]-para[ndp+56]*dat3[,15]-para[ndp+57]*dat3[,16],mean=0,sd=para[ndp+58],log=TRUE))+sum(dnorm(dat3[,11]-para[ndp+59]-para[ndp+60]*dat3[,1]-para[ndp+61]*dat3[,14]-para[ndp+62]*dat3[,15]-para[ndp+63]*dat3[,16],mean=0,sd=para[ndp+64],log=TRUE))+sum(dnorm(dat3[,12]-para[ndp+65]-para[ndp+66]*dat3[,1]-para[ndp+67]*dat3[,14]-para[ndp+68]*dat3[,15]-para[ndp+69]*dat3[,16],mean=0,sd=para[ndp+70],log=TRUE))+sum(dnorm(dat3[,13]-para[ndp+71]-para[ndp+72]*dat3[,1]-para[ndp+73]*dat3[,14]-para[ndp+74]*dat3[,15]-para[ndp+75]*dat3[,16],mean=0,sd=para[ndp+76],log=TRUE))
  
}



out112=fdHess(out111$par,likhess)
H111=out112$Hessian/length(dat3[,1])
(sqrt(diag(inv(-H111)))/sqrt(length(dat3[,1])))